Timothée Leblond
email: timothee.leblond@irt-systemx.fr

Pierre Froment
email: pierre.froment@irt-systemx.fr

Paul De Nazelle
email: paul.denazelle@irt-systemx.fr

Reda Sellakh
email: reda.sellakh@supmeca.fr

Philippe Serré
email: philippe.serre@supmeca.fr

Gaël Chevallier
email: gael.chevallier@univ-fcomte.fr

Gradient-based Optimization of Parameterized CAD Geometries

Keywords: CAD parameters, sensitivity, optimization, CFD, adjoint solvers

This paper presents a novative method for gradient-based optimization with regard to CAD parameters. This method allows to respect manufacturing, design and surface quality rules, particularly required in computational fluid dynamics (CFD). It can be used with any objective function available in adjoint solvers (both in structural and CFD). To prove the efficiency of the method, the workflow schedule was apply on an air-conditionning duct so as to maximize the eigen frequency and then to minimize the pressure drop. 2.

Introduction

Thanks to HPC, numerical optimization methods are more and more used to determine an optimal shape at a lower cost in faster. Among these methods, both design of experiments and surrogate modeling methods allow to work directly on CAD parameters with the exception that the design space, i.e. the number of parameters, should be confined enough in order to be explored within a reasonable computational time. Other methods are based on the gradients provided by adjoint solvers, that is to say on the sensitivity of a cost function with respect to the displacement boundaries. The sensitivity is used to know how to change the shape at any node of the mesh to obtain better performance. They naturally get over the number of degrees of freedom but they are computed with respect to the coordinates of the vertices of the surface mesh. However, since manufacturing constraints are difficult to express mathematically, they are not taken into account during the mesh morphing. This drawback can be avoid by using a CAD model that implicitly includes these manufacturing constraints. To take advantage of both approaches, this method proposes to extend these gradients to CAD parameters. This paper presents the developments of this method and its integration into an optimization loop. Applications on structural and CFD problems will be presented to prove the feasibility of this approach and the possible gains about the number of resolutions.

Shape Sensitivity with regard to CAD Parameters

During an optimization, the shape of the geometry is modified in order to improve its performance regarding an objective function. Here, it can only be modified thanks to CAD parameters. So, it is necessary to determine the shape sensitivity of the objective function with regard to CAD parameters (d α J) given by Eq. [START_REF] Armstrong | Linking adjoint sensitivity maps with CAD parameters, Evolutionary methods for design, optimization and control[END_REF].

d α J = ∂ x J × ∂ α x (1)
∂ x J represents the shape sensitivity of the objective function J with regard to node displacements. It is given by adjoint solvers and indicates how to move the nodes of the mesh to improve the objective function. In that scheme, any objective function available in adjoint solvers can be applied in this equation. The term ∂ α x represents the shape sensitivity with regard to CAD parameters and indicates how the shape is impacted by a CAD parameter perturbation.

State of the Art

The bibliography have pesented three ways to compute the sensitivity of the shape with regard to CAD parameters. The first one is an analytic way that suggests to derivate the CAD surface equations. Yu et al. [START_REF] Yu | CAD-based shape optimisation using adjoint sensitivities[END_REF] present the main advantage of this method, i.e. the accuracy of the results. However, it is necessary to know the explicite definition of the shape in order to derivate these equations. Mostly, these formal expressions are not available in the CAD systems that are often proprietary software and that use implicite geometric solvers. It is therefore difficult to compute shape derivation.

The numerical approach consists in evaluating the gradient by finite differences as proposed by Armstrong et al. [START_REF] Armstrong | Linking adjoint sensitivity maps with CAD parameters, Evolutionary methods for design, optimization and control[END_REF]. It compares the initial and perturbed geometries prior to generate the meshes. This approach requires the disponibility of the geometric representation of CAD models [START_REF] Hardee | A CAD-based design parameterization for shape optimization of elastic solids[END_REF] so as to be able to work with u, v coordinates. Moreover, a parameter variation may result in additionnal patches and different naming [START_REF] Kripac | A mechanism for persistently naming topological entities in history-based parametric solid models[END_REF]. These changes makes surface identification difficult [START_REF] Robinson | Sensitivity-based optimization of parameterized CAD geometries[END_REF].

Another numerical approach is based on the discretization of the CAD geometry. The evaluation of the gradient is done by finite differences between both the meshes of the initial shape and the perturbed shape. The main difficulty is to determine the deformation field that projects one mesh on the other and that limits the distorsion of cells. The deformation of the initial mesh is presented in many methods and Toivanen et al. [START_REF] Toivanen | A new method for creating sparse design velocity fields[END_REF] present one of them but the cell distorsion is hardly mastered. Robinson et al. [START_REF] Robinson | Sensitivity-based optimization of parameterized CAD geometries[END_REF] propose to project the nodes of the initial mesh on the deformed one along the normal but it causes difficulties on the process and requires adjustements.

Proposed Method

When meshes are used to compute the gradient, one of the main problems is the absence of isoconnectivity between both the mesh of the initial shape and the mesh of the geometry after the variation of a parameter. The evaluation of the displacement of any node of the mesh is used to compute the sensitivity. To determine this displacement, the mesh of the initial shape is morphed on the mesh of the perturbed shape. During this operation, it is easy to follow each node and create an isoconnectivity between the mesh of both geometries. The harmonic projection [START_REF] Gu | Computational conformal geometry[END_REF][START_REF] Remacle | High-quality surface remeshing using harmonic maps[END_REF] is used to project 3D geometries on a 2D parametric domain, it corresponds to a reparameterization in 2D space. The latter is considered as a reference space for both shapes.

The theoretical process is detailled with a geometry Ω α with one frontier. At each iteration of the optimization, α 0 represents the set of CAD parameters for the initial shape Ω α 0 , considered as the reference geometry. The perturbed geometry Ω α i is obtained after the variation of i th parameter. The respective meshes of the initial and the perturbed shapes are noted Ωα 0 and Ωα i . Each mesh is projected by a harmonic transformation on a unit disk, respectively noted Dα 0 and Dα i . Given U α 0 and U α i the functions used to obtain respectively these disks and U -1 α 0 and U -1 α i the inverse functions:

U χ : Ωχ ⊂ R 3 -→ Dχ (0, 1) (x, y, z) -→ u χ (x, y, z) v χ (x, y, z) χ ∈ {α 0 , α 1 , . . . , α i , . . . , α p } (2)
The figure 1 gives an illustration of the previous step. For example, a cube, opened on the lower face, is considered as the reference shape Ω α 0 . The variation of the parameter controlling the height of an edge provides the perturbed shape Ω α i . The respective meshes Ωα 0 and Ωα i , done separately, are projected onto the disks Dα 0 and Dα i . The harmonic function is one-to-one, so it is possible to apply U -1 α i on the disk Dα 0 . The nodes of Dα 0 are interpolated with ones of Dα i . So, the initial mesh Ωα 0 is fitted on the perturbed mesh Ωα i . Thus, the new values of the coordinates of Ωα 0 are known. The evaluation of the perturbation is given by the finite differences of the old and new coordinates of each node, as given by Eq.(3).

Ωα 0 Dα 0 Dα i Ωα i U αi U α 0
∂ α i Ω U -1 α i Dα 0 (0, 1) -Ωα 0 δ α i i = 1 . . . p (3)
It is now possible to compute the Eq.( 1). Nevertheless, the presented method requires that the 3D surface and the parametric domain are homeomorphic. In this paper, we handle CAD geometries so has to make them be home-omorphic to a disk Dχ . Indeed, the surface is cut to obtain a homeomorphic surface to a disk. This solution was used to verify the feasibility of the proposed approach in an optimization loop.

Applications in structure and CFD

The optimization loop uses the result of the Eq.( 1) into a gradient descent algorithm. The optimization loop can be used in structural and CFD domains. In direct optimization, even if the computation time is limited, it is possible with the gradient to directly improve the solution, particularly in CFD where the computational cost can be very expensive. This method, breaking away from the usual industrial process, was applied on different academic test cases. Two of them are detailled after and a summary of specific ones is done in the table 1.

As presented by the figure 2, the initial values of the parameters α 0 are given to the CAD software. The outputs are the initial geometry Ω α 0 and every geometries Ω α i corresponding to the variation of each parameter α i . In this case, the variation of each parameters are chosen arbitrarily at 5% of the initial value. All these geometries are meshed and the mesh of the initial one is given directly to the primary and the adjoint solvers. During the resolution, every meshes, corresponding to each perturbed geometries, are compared to the initial mesh by the processus presented in 4.2. Both information given by the adjoint solver and the sensitivity computation is used in an optimizer to determine the new values of α 0 for the next iterations. First, gradients from Nastran (available since 1995) were used to prove the efficiency of the method. The optimization problem was to maximize the eigen frequency of an air-conditionning duct. The CAD geometry (figure 3(a)) is a multi-section surface defined by 7 sections. The parameters and the boundary conditions, applied to the first and the last section, give 13 parameters and a frequency of 300 Hz. After about 40 iterations, the frequency was more than 1000 Hz. So, the optimization improve significantly the shape regarding this objective function only by modifying CAD parameters. Second, the optimization problem was to minimize the pressure drop of the air-conditionning duct. The Star-CCM+ adjoint solver was used to compute the sensitivity of the objective function. The CAD geometry (figure 3(b)) has a S-bend profile defined by a multi-section surface with 7 sections. The parameters of each section and the boundary conditions give a model with 25 parameters and an initial pressure drop of 530 Pa. After 14 iterations, the gain was around 40 percent only by modifying CAD parameters.

Ω α 0 Ω α 1 • • • Ω αn Ωα 0 Ωα 1 • • • Ωαn J, ∂ x J ∂ α x α 0

Conclusion

Current results are obtained by cutting the CAD geometries in order to make the 3D surface be homeomorphic to a disk. However, this requires an extra step to the designer to find a cutting surface valid whatever the value of parameters. Futher developments will avoid this extra step and project multi-boundary surface on the disk.

Results prove the feasibility of the presented method with significant gain regarding the number of resolutions. It would be interesting to use these gradients in more complex algorithms such as conjugate gradient or quasi-Newton or in response surface modeling. In this paper, structural and CFD domains are separately treated. Further developments will allow to apply this loop on multiphysic cases.

Figure 1 :

 1 Figure 1: Projection of the mesh of the initial and perturbed shapes

Figure 2 :

 2 Figure 2: Scheme of the optimization loop

Figure 3 :

 3 Figure 3: CAD geometries used in optimization loops

Table 1 :

 1 Summary of results obtained with the optimization loop in structure and CFD -Structure 1 : maximizing the eigen frequency -Structure 2 : minimizing this displacement due to a force -CFD: minimizing the pressure drop

	Physic	Profile	Number of	Number of	Objective function	Gain
	involved	of duct	parameters	iterations	initial	final	[%]
	Structure 1	Bend	13	36	306.2 Hz	1007.6 Hz	229
	Structure 1	Bend	14	52	217.6 Hz	673.9 Hz	210
	Structure 2	Bend	11	45	9.05 mm	1.21 mm	86.6
	Structure 2	Bend	30	66	5.68 mm	0.55 mm	90.3
	CFD	Bend	10	23	4.86 Pa	2.44 Pa	49.8
	CFD	Bend	25	15	4.32 Pa	1.75 Pa	59.5
	CFD	S-bend	25	14	532 Pa	302 Pa	43.2

th -12 th , June 2015, Sydney Australia

Acknowledgements

This research work has been carried out under the leadership of the Technological Research Institute SystemX, and therefore granted with public funds within the scope of the French Program "Investissements d'Avenir". We thanks CD-Adapco France for providing us the Star-CCM+ licenses.


A Bezerianos

F Chevalier

P Dragicevic

N Elmqvist

J D Fekete

GraphDice: A System for Exploring Multivariate Social Networks

Keywords: (according to ACM CCS): H, 5, 1 [Information Systems]: Information Interfaces and Presentation-User Interfaces; E, 1 [Data]: Data Structures-Graphs and Networks

Figure 1: Exploration of the InfoVis 2004 Contest co-authorship dataset using GraphDice. On the left is the main visualization window of GraphDice including (a) an overview plot matrix, (b) a selection history tool, (c) a selection query window, (d) a main plot, and (e) a toolbar. Overlapping nodes in the main plot are drawn using jitter (visible in the yellow selection query). On the right are actor (f) and link (g) tables with query data entries highlighted in the corresponding color.

Introduction

A quite recent development in social network analysis (SNA) [START_REF] Wasserman | Social Network Analysis[END_REF] has been the adoption of visualization to explore networks and support social scientists in detecting, understanding, and characterizing unexpected patterns and trends in complex social networks [Ada06, HF06, HF07]. real social networks collected and built by sociologists and historians are rich in attributes for both actors (graph nodes) and relations between actors (graph edges). To compound the problem, the standard network analysis algorithms for computing degree, centrality, and graph layout that these domain experts routinely use create additional node and edge attributes. The process of analyzing such networks is called multivariate social network analysis.

To address the dearth of multivariate network visualization tools, we present GraphDice, a tool based on the Scat-terDice [START_REF] Elmqvist | Rolling the dice: Multidimensional visual exploration using scatterplot matrix navigation[END_REF] dimensional navigation method, suitably adapted for network visualization. GraphDice (Fig. 1) is the first tool to use a plot matrix to navigate multivariate graphs. It is designed to be simple to learn and use by social science researchers, while providing very rich navigation and exploration capabilities. Further contributions of this article include the unified visualizations of edge and node attributes, by supporting interval values on the nodes and projecting link attributes on the nodes as intervals (Fig. 2), as well as interesting visual configurations provided through computed attributes to unify the navigation.

In this paper we first discuss previous work on multivariate social network analysis and visualization. We then describe a set of design goals followed during the design of GraphDice, followed by an overview of the system. Finally, we present a number of real datasets visualized by our tool and discuss results from a user feedback session.

Background

In this section we review work on social network and general multidimensional visualization, including for graphs.

Social Network Analysis (SNA)

A social network is a graph structure whose vertices are social entities such as persons or organizations, and edges are relations between these entities. Social network analysis (SNA) [START_REF] Wasserman | Social Network Analysis[END_REF], correspondingly, consists of finding structural properties in the network, mostly through statistical and graph-theoretical methods. Social networks typically associate data to vertices and edges. For example, if vertices are persons, they can have a name, a birth date, a position in a company, and many more attributes. Similarly, relations can also have attributes, such as date of friendship connection.

Just like statistics, SNA can be confirmatory or exploratory. Most books on SNA are devoted to confirmatory methods, where there is an initial hypothesis about the structure of the network, and several methods are described for how to confirm this hypothesis. However, exploratory SNA is becoming an increasing opportunity due to the availability of vast amounts of data, either from online social networks, like Facebook, or constructed from social sites, such as Wikipedia or Flickr. Moreover, organizations increasingly maintain information about their members in databases that are easy to explore as social networks by focusing on any relation such as the organizational chart, project history, or email exchanges. These social networks are very rich, and besides confirming a-priori hypotheses, exploring these networks often reveals interesting and unexpected phenomena and structures. This ready availability of data also means that more people with no formal training in social network analysis have access to networks which are meaningful to them, such as relationships between cinema actors, relationships between their family and friends. However, most established tools for analyzing social networks, such as UCINET [START_REF] Borgatti | UCINET 6 for Windows: Software for social network analysis[END_REF], are intended for confirmatory analysis, and are therefore weak at exploration and require substantial training. The International Network for Social Network Analysis [INS] lists about 60 software tools that fall into three categories: confirmatory analysis, exploratory analysis and network visualization. The popular Pajek graph visualization system can be used for exploration [dNMB05], but not for interactive information visualization because each exploration step is very indirect, involving long menus, several sub-menus, or large dialog boxes, once again requiring special training.

Computed attributes play an important role in SNA, like centrality or roles (e.g. hub, gate-keeper). Graph-theoretical properties, such as degree of a vertex (in-degree and outdegree for directed graphs) or clustering coefficients, also provide important information about the status of actors in the network. However, most measures have sometimes numerous variants, (e.g. centrality has tens of variants that, depending on the network structure, represent more or less faithfully the intuitive notion of centrality). Therefore, even when analysts want a measure of centrality, they often compute many of them and check -statistically or visuallywhich best fits their intuition about the reality of the world.

In practice, standard SNA consists of computing useful measures on the network and finding structures based on the topology, the intrinsic attributes and these measures. This is why most analysts -whether professional or casual -often need to explore multivariate networks with typically 10 to 30 measures per actor or relation.

Tasks

There has recently been considerable work on enumerating important tasks for social network analysis and exploration. Lee et al. [LPP * 06] present a general task taxonomy of: low-level tasks, topology-based tasks, attribute-based tasks, browsing tasks, and overview tasks.

For social networks, Aris and Shneiderman [START_REF] Shneiderman | Network visualization by semantic substrates[END_REF]AS07] propose a set of 6 challenges and 10 tasks. The challenges are about visualizing: C1) Basic networks with nodes and links; C2) node labels; C3) link labels; C4) directed networks; C5) node attributes; and C6) link attributes. The tasks listed start from low-level and add more SNA related tasks: T3) for every node, find the nodes that are distance 1, 2, 3 . . . away; T4) for every node, find betweenness centrality; T5) for every node, find structural prestige; T6) find diameter of the network; T7) identify strongly connected or compact clusters; T8) for a given pair of nodes, find shortest path between them; T9) for every node/link, read the label; and T10) find all nodes/links with a given label/attribute. Finally, Henry et al. [HFM07] list three high-level tasks specific to SNA: U1) identify communities; U2) identify central actors; and U3) analyze roles and positions.

Social Network Visualization

Social networks are usually visualized using node-link diagrams or adjacency matrices. Several social network systems have been designed, including VisOne [START_REF] Brandes | Graph Drawing Software[END_REF], Vizster [HB05] and MatrixExplorer [HF06].

Most general graph visualization tools listed can be used to visually explore social networks but are not specialized for this purpose. They provide graph-layout algorithms and the possibility of assigning data attributes to visual attributes such as color and size, up to three or four attributes, in addition to positions. However, they may lack specialized layout methods for social networks, such as [Noa05].

Systems targeted towards SNA compute a layout based on the topology -sometimes taking some attribute into account in the layout computation -and allow the use of other attributes to decorate the nodes and links. This effectively adds 3 to 5 extra visualized attributes before the visualization becomes too complex (width, height, background color, outline color, transparency). Since social networks typically have 10 to 30 attributes associated with the nodes, exploring their values becomes tedious.

Furthermore, these visualization tools focus on the topology of the networks, not on their attributes (C5,C6). Topology is clearly very important (C1,C2), but for social networks, other perspectives may be as important to discover correlations between attribute values. Therefore, some recent systems have tried to offer better support for attributebased navigation in multivariate networks.

Multivariate Network Visualization

Since 2006, there has been an increasing interest in visualizing multivariate networks in more clever ways than changing visual attributes of nodes and links. Wattenberg's Piv-otGraph system [Wat06] relies on the OLAP (Online Analytical Processing) database model to aggregate network vertices and drill up or down according to categorical attribute values. Initially, the PivotGraph system shows one node that represents the whole dataset aggregated. The user can drill down using any categorical attribute, causing the initial point to be split into several points -one for each value of the attribute -that are spread evenly on the horizontal axis. The point sizes are proportional to the number of individuals (database rows) holding the specified value; links are drawn between the points with a width proportional to the number of edges connecting the points. A second attribute can be used for splitting the vertical dimension; the points are then laid out in a regular 2D grid. PivotGraph can visualize very large databases along any categorical dimension as long as the number of categories remains small enough (C5,C6,T9,T10). However, it is not meant to address the tasks and challenges of full-scale SNA.

Following-up on the idea of exploring multiple attributes of a network, semantic substrates [SA06, AS07] use a different approach with no aggregation but several panes on the same window. The technique displays a social network showing several attributes unfolded. Each pane visualizes two attributes of the vertices using two dimensions. Various dynamic query widgets are placed on the window to allow filtering. Furthermore, edge filtering can be done on the panes by drawing a box around the vertices of interest, hiding links outside the box. While substrates support most of the challenges C1-C6, building the panes requires an initial setup phase. This phase needs computer science skills beyond that of social science researchers and, even more difficult, a clear understanding of the important attributes to show and how to lay them out in space to reveal important relationships.

Pretorius et al. [Pre08, PvW08

] describe several methods to visualize multivariate networks, mainly applied to transition networks: brushing and linking of a node-link diagram with parallel histograms for each dimension, linear graph visualization of a hierarchical clustering of vertices with various attributes superimposed, and hybrid parallel coordinates. All these methods seem promising but, except the first, they are also quite complex, both visually and cognitively. The list of tasks supported by these visualizations do not focus on structural groups (T7,T8) but on attributes; thus, it remains to be seen whether they can be used effectively for SNA tasks.

In summary, visualizing multivariate networks is recognized as an important but difficult challenge. PivotGraph remains simple but only shows categorical attributes and does not support most SNA tasks. Semantic substrates support the identified tasks, but are complex to set up and only support visualization of a small set of attributes at a time. Pretorius' visualizations are complex, difficult to learn, and not meant to support the group tasks required by SNA. Our challenge is to provide a simple visualization tool supporting a more complete set of the SNA tasks.

Design Rationale

Henry and Fekete [HF06] list 13 requirements for SNA, as collected from participatory design sessions. We take into account overall system design requirements: R3) show an overview; R6) use analytical information (SNA measures); R7) preference for interaction (direct manipulation) vs. parameter tuning; R8) layout automatically and interactively; c 2017 The Author(s) R9) filtering is necessary. Particularly important for this article are requirements related to multiple attributes: R5, R12 -Attributes & Consensus: Taking attributes into account makes the difference between graph drawing and information visualization. Participants were not interested in displaying a unique graph; they wanted to build several representations according to the different attributes of the edges and vertices (R5). They also wanted to compare and identify a consensus among actor clusters across representations (R12). The structure of the graph may differ depending on the chosen attributes. Comparing these structures, understanding why they are similar or how they differ was a major concern. Users needed a visualization system which helps them to choose visual variables for each attribute and create multiple views of their dataset. Consulting details for each vertex or edge was also a primary interest and therefore, details should always be visible or quickly accessible.

Multivariate social networks are complex structures and the tools available for their exploration are themselves complex and require substantial training to learn and understand their features, making multivariate SNA a very challenging task. In GraphDice we thus attempt to efficiently support R5 & R12. Our goal is to simplify exploration, feature understanding, and training, by providing a tool that is simple to learn and use, yet powerful enough to analyze networks with hundreds of nodes of up to 30 dimensions. To achieve this goal, we rely on a set of design principles derived from the above requirements, as well as ones that can be seen as a stronger interpretation of the well known rules of direct manipulation [Shn97] and dynamic querying [START_REF] Shneiderman | Dynamic queries for visual information seeking[END_REF]: P1 Use direct manipulation (R7). Allow users to interact directly on the visualizations instead of using complex widgets that take time to understand, and act indirectly on items of interest; limit the number of hidden functions only reachable through menu navigation that disrupts the user's main focus. P2 Keep the widgets simple and small in number. Devote most of the screen to the data, and not on widgets that consume screen real-estate. P3 Keep the number of visual representations small to shorten learning time and limit cognitive fatigue; P4 Keep the interactions and visualizations consistent, to facilitate learning and invite exploration; P5 Visualize the navigation and interaction history (undo/redo states) to facilitate learning and invite traceable, reversable and thus inexpensive exploration. P6 Support (R6) the computation of basic SNA metrics (e.g. centrality, degree, etc). P7 Provide overview, details and attribute filtering (R3,R5). P8 Maintain visual groupings across attributes (R12).

The GraphDice System

GraphDice is a visualization system for exploring multivariate social networks, following the design guidelines of Sec. 3. It extends the ScatterDice [START_REF] Elmqvist | Rolling the dice: Multidimensional visual exploration using scatterplot matrix navigation[END_REF] tool, that was originally designed for navigating and exploring multidimensional tables visualized as multiple scatterplots. Nodelink diagrams can be seen as extensions of scatterplots where points are connected with links. However, the application of scatterplots to network data implies an increase of complexity of its data model that requires different features and design choices than for 2D scatterplots.

The GraphDice system (Fig. 1) consists of a visualization window and data table views. The visualization window in turn includes several components discussed next. Images come from real datasets described in detail in Sec. 6.

Overview Plot Matrix

Like ScatterDice, the overview plot matrix displays a crossdimensional node-link plot for every combination of attributes and arranges them in a large plot matrix [START_REF] Fisherkeller | PRIM-9: An interactive multi-dimensional data display and analysis system[END_REF]. All attributes are assigned a row and a column and a plot thumbnail is drawn at their intersection (Fig. 1a). Scatterplot matrices have never before been used to visualize graphs.

The main plot can be changed by navigating in the plot matrix, similarly to ScatterDice. This navigation is restricted to orthogonal movement along the same row or axis in the matrix: one axis is preserved while the other changes using an animated transition. Instead of simply interpolating the position of each data point, a 3D rotation is performed, as if rotating a cube or rolling dice. These transitions allow the user to extract structure from motion [Ull79] and compare the two dimensions. It also reinforces the metaphor of navigating in the attribute space defined by the plot matrix (P4).

As an extension to ScatterDice, user selections have been made visible on the thumbnails: each selected node is displayed with its selection color. This feature is useful for tracking selections across dimension pairs, and finding attributes where selected nodes remain grouped (P8).

Main Plot

Node-link diagrams can be seen as extensions of scatterplots where data points are connected with links. Therefore, the use of node-link diagrams is consistent with the ScatterDice "dice rolling" paradigm. Because GraphDice allows to focus on a single visual representation at a time (Fig. 1d), the system minimizes learning time and cognitive fatigue (P3).

Attribute Drawing

In multivariate data, attributes are often categorical (e.g. Gender or Country), resulting in multiple nodes overlapping on one or more dimension axes. GraphDice uses jitter [AS07] to provide a clearer visual indication that multiple nodes occupy the same space in the graph. Overlapping nodes are placed around a circle (Fig. 1d), centered on their original location, with a radius proportional to the number of overlapping nodes. Using this jitter mechanism, users can identify at a glance dense data regions. Moreover, they can see the connections between these overlapping nodes as the circular layout provides space to draw the interconnection links. When showing categorical axes, the resulting plot is very similar to a PivotGraph [Wat06] but the aggregation is only visual, so all details remain visible.

Axes can also encode intervals, i.e. value ranges. Intervals allow aggregatation of numerical attributes of actors, e.g., the years of a researcher's publications can be replaced with the year of her first and last publications. GraphDice visualizes intervals as straight lines with short hooks in the ends, to distinguish them from links between actors (Fig. 2).

For edges with a numerical or an ordered categorical attribute, GraphDice automatically creates an associated interval node attribute and aggregates the edge values at each node. In the co-authoring network of Fig. 2, edges are articles with a publication year. GraphDice automatically creates an interval node attribute for the "year" attribute named "e:year". When this attribute is set as an axis in the plot, the links are drawn starting and ending at a position computed from their attribute value. So if "e:year" is assigned to the horizontal x-axis, authors are drawn as straight horizontal lines showing the intervals of their publication years. Links are drawn as curves with their endpoints positioned according to the value of their "year" attribute. Thus GraphDice allows the exploration of node and edge attributes in a unified way (P3, P4). To benefit from this feature, GraphDice supports multiple (parallel) edges between the same nodes.

Link Drawing

Links are drawn using splines rather than straight lines, as it has been shown to improve the legibility of dense nodelink diagrams [FWDP03, AS07]. Directed links are drawn using biased splines as in [FWDP03], where bias and orientation express the link direction (Fig. 3). For undirected networks, GraphDice attempts to find the best link drawing method using a set of simple heuristics for each configuration in the plot matrix. In diagonal scatterplots and 1-D Figure 3: Link directionality in a migration network of cities in 2 discrete administrative departments. The excentric label is used to focus on a single node and its links.

scatterplots links are drawn as semicircles, like in Arc Diagrams [Wat02]. When one of the two axes is categorical, ordered (e.g., rank, x, and y), or is an interval, links are splines bent towards the orthogonal direction (Fig. 2). When both axes are categorical, a mix of vertically and horizontallyoriented splines is used to minimize link overlapping.

The user can further tweak the rendering of individual networks in the matrix by changing the type, direction and amount of link curvature. Finally, links are colored by drawing a gradient between the two endpoint node colors.

Link Animation

To aid navigation between different attributes and help users maintain a mental model of the dataset, every plot maintains the links between nodes, even during plot transitions, irrespective of the visualized attributes. To our knowledge, most multivariate visualization systems for social networks do not provide persistent link presence between attribute views.

As in ScatterDice, GraphDice uses 3D rotations to animate between visualizations: a 2D network is extruded in 3D, rotated, and flattened back to a new 2D network. In Fig. 5c the plot is in a partial transition stage between two attributes. Nodes are animated the same way as datapoints in ScatterDice. In the case of interval attributes, the animation involves a progressive expansion/shrink of single points (nodes) to straight lines (intervals).

Animating links is technically more challenging than nodes, as it is not possible to reconstruct arbitrary links in 3D. E.g., there is no 3D curve that appears as an upward semicircle in a given projection and as a downward semicircle after a 90 • rotation about the vertical axis. To address this issue, we devised the following technique: before the animation starts, we compute consensus 3D curves from the initial and final 2D links (Fig. 4). During the initial plot extrusion phase, 2D links are morphed to consensus 3D curves by linearly interpolating the curves. At the rotation phase, the consensus 3D curves remain unchanged. Finally, at the flattening stage, 3D curves are morphed into the new 2D links.

GraphDice employs two different methods for computing the consensus 3D curves. When the initial and final 2D links have the same direction of curvature (e.g., upwards), they are both extruded in 3D; the two parametric curves obtained are then averaged to form the consensus 3D curve (Fig. 4a). Since the initial and final links have a similar shape, the visual transition to and from the 3D curve appears natural. When the two extruded curves are confounded, this technique reduces to a regular 3D curve reconstruction (Fig. 4b).

In case the initial and final 2D links have different directions of curvature, the above technique produces complexlooking 3D curves with inflection points. We found that using straight lines as consensus 3D curves produces animations that are easier to follow. Therefore, when initial and final curvatures do not match, the link is straightened during the extrusion phase, after which a 3D line rotation is performed, and finally the line bends again as the cube flattens.

Computed Attributes

Before loading a network, we compute additional attributes for each vertex that are useful in network exploration (P6). These attributes are either structural, layouts or orders.

Structural attributes include the node degree (or in and out for directed networks) and betweenness centrality [Fre77] and are computed using the JUNG network library [START_REF] Fisher | JUNG: Java universal network/graph framework[END_REF].

For the layout, we assign the X and Y attributes with positions computed using the edge-repulsion LinLog graph layout [Noa05] (default view). When the network already contains an x and y attribute (lower case), we use it as the initial layout before applying the force-directed layout algorithm. Therefore, networks containing geographical coordinates have a layout as consistent as possible with the expected geographical positions.

Orders are attributes of consecutive values, from 1 to the number of nodes in the network. When the attribute is assigned to an axis, all the nodes take an evenly distributed position along that axis following that particular order. Orders are useful to visually analyze the network. We compute an order that minimizes the bandwidth of the network using the reverse Cuthill-McKee algorithm [CM69]. This ordering groups connected nodes since the link lengths are kept short, and can augment any attribute with a meaningful topological order. Since GraphDice creates an interval node attribute for each edge attribute, we create an order that shows a meaningful progression of the intervals, sorted by increasing average values first, then by increasing range length (Fig. 2).

Contrary to ScatterDice, the order of columns and rows in the overview plot matrix is meaningful and computed automatically as a function of the similarity between dimensions [START_REF] Ankerst | Similarity clustering of dimensions for an enhanced visualization of multidimensional data[END_REF] (except for the X, Y layout dimensions that are placed at the end of the plot matrix). This feature was a request of our social scientist user group who are always interested in the correlations of their attributes. The attribute order can be manually changed through interactive drag-and-drop of rows and columns. Details on how to navigate within the plot matrix by selecting or dragging on plot thumbnails can be found in the ScatterDice article [START_REF] Elmqvist | Rolling the dice: Multidimensional visual exploration using scatterplot matrix navigation[END_REF].

As this paper focuses on visualizing rather than computing network attributes, our prototype automatically computes only a small subset of the structural and layout algorithms of JUNG [START_REF] Fisher | JUNG: Java universal network/graph framework[END_REF]. But as it is important for the user to be able to choose among multiple layout and structural algorithms, we plan to fully integrate JUNG in GraphDice.

Visual Queries on Edges and Nodes

Beyond basic navigation functionality, the original Scatter-Dice tool also provides advanced visual queries/selections on nodes (P1), called query sculpting, using bounding volumes that allow the user to iteratively filter the dataset. The process consists of simple lasso-selecting data items in the main scatterplot using 2D bounding volumes (boxes or convex hulls) and then iteratively refining the selection from other viewpoints while navigating the plot matrix (similarly to high-dimensional brushing [MW95]). The active query hull is blinking to highlight nodes affected by the refinements. Several color-coded queries are available to the user, visible in a query widget (Fig. 1c), much like a Photoshop layer pane. Queries are sculpted in the data space of the current dataset, two dimensions at a time (the dimensions currently viewed in the main plot) and the selection may either be a union or an intersection with the existing items in the current query. The visual representation of queries conforms to the overall navigation metaphor and is thus also "rolled" as a 3D convex hull when the user navigates in the scatterplot matrix. This helps the user to perceive correlations between the extents of a query in two adjacent plots (P12). Queries behave in the same way on both edge and actor axes, providing a consistent direct manipulation approach (P1, P4). Nodes are drawn as small points with a color related to their selection in a visual query (described below) and an intensity mapped to a user selected attribute (initially constant). When a node is selected, it and its edges belong to one or several query layers -each with its unique color. When performing a visual querying in an edge interval attribute view, only edges that belong to this selection are part of the query, making it easy to focus on patterns of a specific edge value, such as a time period. Each query layer shows the proportion of the nodes and links over the number of all actors and links in the network. Contrary to ScatterDice, we color the nodes according to the layers they belong to; accordingly blended if a node belongs to several layers. A selection layer is visualized as the convex hull of all the nodes it contains, colored with the layer's color 90% translucent.

Query History

Maintaining the history of the exploration process (P5) is important to help users revisit important views and keep track of their exploration [HMSA08]. In the main window, we have added an interactive query history tool (Fig. 1b) to help users in revisiting previous steps of exploration. Each time a new query is formulated by the use of the lasso or directly in the table, the current state of the selection and the current axes assignment are stored as a new query step in the history panel. The thumbnail of the current scatterplot is added to the history panel as being the last (and current) query state, showing the current state of the selection over the thumbnail. The thumbnails with the same selection as the current ones are outlined. If the new selection corresponds to the previous recorded step, we consider that the new query corresponds to an undo action: we remove the last step from the history and revert to the previous step; if the current selection state is not the last in the history when a new query is formulated, the states after the current one are removed.

Each thumbnail is a three-state button. While hovering over a thumbnail, the selections in that history state are shown transiently on the main plot. Clicking on a thumbnail permanently installs the selections on the current attribute plot. Clicking again on that thumbnail "rolls the dice" to its associated stored plot view. The history panel thus allows for going back to a previous state, both in reinstalling the selections in a first step, and then going back to the attribute plot view where the queries were performed.

Other functionality

As in ScatterDice, users can navigate in the overview plot matrix by clicking on individual cells or by pressing the arrows keys. In GraphDice, users can additionally navigate from within the plot view by invoking a control menu [PLVB00] (Fig. 5). This menu allows them to select both the axis to "roll" (either x or y) and the new dimension to display. Releasing the mouse button on a menu item triggers the animation, and dragging past a menu item pop-ups a control slider that allows precise control of the 3D transition.

A toolbar is also provided for turning on and off different visualization options (Fig. 1e). Options include auto-zoom, labels, and the display of in-or out-links within selections. The toolbar supports fast pre-visualization by using threestate buttons similar to the 3 state history buttons. Finally, an excentric labels lens [FP99] shows details of specific nodes.

Data Tables

As a supplement to the main visualization window showing node attributes, we provide users with two tabular views (Fig. 1 right) of their dataset (nodes and attributes/dimensions) and a table of their links. Interacting with the dataset table is similar to interacting with a spreadsheet: users can sort data by dimension, select multiple entries, rearrange dimensions, perform text searches, etc. By default, this table includes all data, but as data gets selected by users on the main visualization, the table entries are updated to reflect or adjust the current user selection. More precisely, when the excentric label is active, the table only shows the nodes inside the lens. When zoom is active, the table only shows the selected nodes.

Thus, the tables act as a details-on demand mechanism by providing information on entries (P7), as well as a verification and correction mechanism for user data selections. Finally, the data table is used for dynamic queries: when rows are selected, a right click triggers a popup menu to choose a query layer, assigning the selected rows to this layer.

Following Guidelines and Supporting Tasks

We designed GraphDice following the principles of Sec. 3, as indicated in our system description. Although our system supports all the described principles, we kept the number of widgets small and consistent in functionality, concentrating more on having space for the main network plot (P1, P2).

Our system is designed to address a large number of challenges and tasks described in Sec. 3. We support visualization of the network (C1) , and node and link label visualization (C2,C3) through our excentric label, 3 state buttons and interactive data tables. By using curved links we display directionality (C4), while our plot matrix visualization and plot transitions allows us to display a large number of node and edge attributes (C5,C6). Moreover, our data table is used to implement the requirements C5 (show node attributes).

When it comes to tasks, we support betweenness and centrality finding (T4, U2) through our computed attributes and relevant plots. Moreover, users can visually identify strongly connected and compact actor clusters based on proximity in our layout algorithms (T7,U1), and compare them across dimensions using visual queries and plot transitions. By using either excentric label or visual queries, node/link labels are readable (T9). Finally, our data tables act as a details-on demand mechanism (providing information on entries) for task T10 (find all nodes with a given attribute), as well as a verification and correction mechanism for user data selections.

Datasets

We built GraphDice to aid visual exploration of multivariate social networks. To evaluate its utility, we used it to visualize 3 real datasets, that differ in number of nodes, edges, node dimensions and connection patterns. Our goal was to determine if dataset characteristics and patterns were visible using GraphDice as an exploratory visualization tool. The images referenced show examples of dataset views and the mentioned exploratory findings. While our paper does not focus on scalability, we have successfully rendered in realtime graphs of up to 2000 nodes and 6000 edges. Visualizing larger graphs while allowing readability would require aggregation, an issue that is out of the scope of this paper.

InfoVis Co-authorship Dataset

The first exploration dataset is the largest connected component of the InfoVis coauthorship dataset from 1998-2002, extracted from the InfoVis 2004 Contest [FGP04]. It contains 694 publications (edges), among 149 authors (nodes), categorized in 10 dimensions. Of these dimensions, 5 are the main dataset attributes (papers, citations, first publication date, last publication date), while the rest are computed automatically by the system (centrality, rank, degree, constant, as well as X, Y coordinates for a 2D graph layout). One edge dimension (publication year) was present.

We focus on part of the network (Fig. 5a), by zooming into query selections representing 2 closely connected groups of co-authors (blue, green) and their connections (red). We observe through centrality (Fig. 5d) that the first group has a single central actor, indicating student-advisor relations, whereas the second has multiple such actors pointing to a more connected research group, like a research organization. The 2 different groups are linked to each other only through connections between these major actors.

Organization Collaboration Dataset

The second dataset (Fig. 2b) is the collaboration network in terms of document creation, between the research teams of a large nation-wide research organization. The dataset consists of 197 teams (nodes) and 1581 collaboration articles between teams (edges). There are 8 dimensions for the teams: system computed ones (as above) as well as creation year and publication count. There are 2 dimensions for the publications (edges), year of publication and count (how many such publications between two teams within a year).

By observing the co-authorship trends of different teams across multiple years, we see an increased amount of publications in 2006 (Fig. 2b). What seems to differentiate this year from others is that a large number of teams disbanded because their project finished (end of interval associated with teams in the y-axis). Focusing on these projects (Fig. 2a), we see that the number of publications the in year they disbanded increased from previous years. We hypothesize that as this was the last year of the projects, a large number of research conducted in previous years was finally published.

Historical Migration Dataset

The third dataset comes from a historian specializing in SNA of migration patterns in 19th century France. It consists of 1776 population migration paths (edges) between 75 cities (nodes) in northern France in 3 time periods. The 18 dimensions of the dataset included, but were not limited to, dominant language of city (French or Flemish), percentage of population occupation (worker, farmer, etc), proximity to riverbanks, administration department, longitude and latitude, etc., as well as automatically computed measures (as above). An edge attribute (migration year) was also present.

By observing the migration network using GraphDice in our user feedback session, the historian was able to form a new hypothesis (Fig. 3). Migrations to and from Godewaersv and cities within the same department are balanced and more frequent, but mainly outgoing migrations take place between this city and cities in the other department.

Preliminary User Feedback

We conducted a full day workshop with a historian specializing in SNA. Although a more formal user study is required, this session helped us gather informal preliminary user feedback and examine wether GraphDice can be used for realistic SNA tasks. The participant used GraphDice to explore her historical migration dataset (Sec. 6.3) and was able to use the tool effectively with less than 15 min. of training.

Our user immediately commented on the usefulness of the plot matrix with the miniature plots, stating that it is useful to "have all combinations of attributes immediately available to quickly explore hypotheses". This was more pronounced when queries were active as "the query feedback in the overview matrix reveals interesting patterns and then guides the exploration as it stirs up unexpected questions.". She found the default X,Y view very useful as "it allows to quickly get back to the traditional representation of the network to rebuild my referring mental map". She also mentioned that having the matrix always visible made it "easier to remember which combinations she had already viewed". We thus feel we have provided users with a successful representation of the dimension space that is easy to navigate. Animated transitions between dimensions were deemed useful as a means to "compare similar attributes". By interactively controlling the transition between two attributes (e.g. two centrality measures) she could "quickly dismiss one from further analysis because they are too similar".

Queries were deemed very useful in "quickly exploring correlations between attributes". Using the convex hulls of a selection in one dimension, the historian was able to explain the selection by attributes in other dimensions. She also commented that suppressing links inside and outside a selection helped in the understanding of inter-community patterns. The synchronized data tables were also deemed extremely useful, as they "provide details and security and confidence that you are looking at the correct data".

As for the drawing of biased splines to show direction, she enjoyed "being able to see immediately reciprocity of migrations between pairs or groups of nodes", something she was "not able to see with the classical arrow representation".

Through the edge interval view tied to different time periods, she was able to see the evolution of the different cities across time. She found it very easy to identify stable and changing patterns between the three periods, as they are shown concurrently without interfering with each other. She added that jitter helped to quickly identify clusters and decide which attribute values to focus on (e.g. medium cities at river-banks), while the animated transition between dimensions gave more information on the jittered nodes. She also requested additional features for graph annotation. She found the query history useful for re-tracing exploration steps, but would also like to take snapshots of "preferred moments" (including the cursor position) and store them in a history of favorites, so as to quickly revisit views of interest and replay scenarios. Being able to mark preferred views in the matrix was also mentioned: she visually marked one graph by customizing link curvature to make it distinct, but a more light-weight mechanism was requested.

Our historian commented on how easy it is to manipulate GraphDice compared to her other tools, like UCINET or Pajek. She pointed out it is very well suited for discovering patterns and learning new datasets, as it provides visual representations of a network from many points of view. She added that even if she doesn't want to learn statistics (a requirement in most current tools), she can still get results with GraphDice. Lastly she mentioned she would use GraphDice to communicate her findings, as she surmised that a live demo of a well-driven scenario on an easy to use tool, would make it easier to understand and correlate her findings than showing a set of well-chosen snapshots of the datasets.

Conclusions and Future Work

We have presented GraphDice, a tool for multivariate visualization of social networks. The GraphDice tool is based on a comprehensive approach to visualizing attributes -derived and intrinsic alike -using an attribute plot for information display and an overview plot matrix for navigation, similar to the ScatterDice system [START_REF] Elmqvist | Rolling the dice: Multidimensional visual exploration using scatterplot matrix navigation[END_REF]. Attribute plots draw nodes according to their values in the two attributes displayed by the plot, and also show edges between the nodes as visual links. Of course, point occlusion is a major issue when aggregating nodes according to a potentially small set of attribute values, and we address this by introducing node jitter proportional to the amount of overlap in the visualization. We have also showcased usage of the tool through an initial participatory design workshop involving a historian.

Our future work entails extending the GraphDice system to sophisticated history storage and navigation beyond simple undo mechanisms and query storage. We plan to augment our system to allow users to take snapshots of their navigation history as landmarks in their dataset navigation, and to provide notification mechanisms to inform users when they have reached a navigation state that is similar to previously visited states. To further support navigation and long term exploration within the attribute space, we also plan to create annotation mechanisms (such as bookmarking interesting plots in the plot matrix or annotating the main plot). Furthermore, we intend to allow users to compute additional graph metrics, including different clustering and centrality coefficients. We also plan to perform a longitudinal user study of social scientists analyzing their networks using GraphDice.

An online WebStart-enabled version of the system can be accessed at: http://www.aviz.fr/graphdice/

Figure 2 :

 2 Figure 2: An edge e:year attribute (x-axis) projected as an interval node attribute (straight lines along x-axis, 1 per author). The author nodes (y-axis) are ordered so that similar intervals are close to each other. In this co-authorship network (a), we use an excentric label to focus on part of the network (b). Notice hooks at ends of intervals.

Figure 4 :

 4 Figure 4: (a) A consensus 3D curve is built by extruding the initial and final links and averaging them. (b) When links are splines whose intermediary control points are kept unchanged in terms of height and relative x-position, this method reduces to a 3D curve reconstruction.

Figure 5 :

 5 Figure 5: (a) X,Y layout for a graph and queries. Menu navigation: (a) Using the right mouse button the user invokes the dimension menu, with the currently selected dimensions marked in red. (b) As the user drags over a menu the new vertical dimension is highlighted in green. If the user releases the menu, a default speed animation transitions to the new centrality dimension (d). (c) If the user continues dragging, a fast slider appears allowing the user to control the animation speed.

c 2017 The Author(s)

Acknowledgments

We would like to thank Claire Lemercier for her feedback. c 2017 The Author(s)


Pavel Svehla

Khieu-Van Nguyen

Jing-Rebecca Li

Luisa Ciobanu
email: luisa.ciobanu@cea.fr

Quantitative DLA-based compressed sensing for T1-weighted acquisitions

Keywords:

Introduction

Recent advances in the static magnetic field strength of magnetic resonance scanners and in the radio-frequency (RF) detector designs has allowed magnetic resonance microscopy (MRM) to reach spatial resolutions suitable for functional imaging of single cells [START_REF] Jelescu | Highlighting manganese dynamics in the nervous system of Aplysia californica using MEMRI at ultra-high field[END_REF][START_REF] Radecki | Functional magnetic resonance microscopy at single-cell resolution in Aplysia californica[END_REF]. However, in order to reach the full potential of MRM it is necessary to reduce the currently long acquisition times required for obtaining high resolution images. Based on the fact that MR images, among other types of images, are compressible, an image can be reconstructed from a small number of random measurements [START_REF] Donoho | Compressed sensing[END_REF]. This finding opened the field of Compressed Sensing (CS) which can significantly reduce the MRI scan time and found numerous applications in preclinical (4) and clinical (5) imaging.

In CS, high-quality images can be obtained from data sampled well below the Nyquist rate provided that the sampling pattern is incoherent, the images are sparse in a transform domain, and a sparsity-promoting iterative reconstruction is used [START_REF] Donoho | Compressed sensing[END_REF]. The CS method has been previously utilized for the acceleration of T 1 weighted acquisitions for knee cartilage quantification [START_REF] Pandit | Accelerated T1 acquisition for knee cartilage quantification using compressed sensing and data-driven parallel imaging: A feasibility study[END_REF] as well as for MEMRI [START_REF] Li | Fast cardiac T1 mapping in mice using a model-based compressed sensing method[END_REF]. In case of the MEMRI study, CS with random k-space undersampling patterns was employed for fast cardiac T 1 mapping in mice [START_REF] Li | Fast cardiac T1 mapping in mice using a model-based compressed sensing method[END_REF], demonstrating the feasibility and performance of this approach. Both studies used random undersampling schemes in the high frequency domain while fully sampling the low frequency domain, which has been shown to reach a similar performance to that of the polynomial undersampling algorithms [START_REF] Li | Fast cardiac T1 mapping in mice using a model-based compressed sensing method[END_REF]. Recently we have introduced a different approach for undersampling the k-space based on using the diffusion limited aggregation (DLA) random growth model to obtain reduced acquisition patterns in the phase encoding directions. We have demonstrated that this DLA-CS algorithm performs better than the polynomial approach and validated its use in high resolution T 2 weighted imaging. In the present study we implement DLA-CS for T 1 weighted acquisitions in order to perform highresolution quantitative functional MEMRI and we evaluate its performance.

Methods

Undersampling pattern generation

Undersampling pattern generation for DLA-CS RARE (Rapid Acquisition with Relaxation Enhancement) acquisitions was previously reported by Nguyen et al. (8). Briefly, the two phase encoding directions in a Cartesian 3D trajectory were undersampled using an acquisition pattern based on the diffusion limited aggregation random growth model [START_REF] Witten | Diffusion-Limited Aggregation, a Kinetic Critical Phenomenon[END_REF] with the k-space points in the resulting patterns always being restricted to be a subset of the fully sampled Cartesian k-space points.

In this study, following the same procedure, acquisition patterns were generated for seven undersampling ratios ranging from 30% to 90% for a T 1 weighted FLASH (Fast Low Angle Shot) acquisition. For each undersampling ratio, 300 sets of undersampling patterns, each consisting of 100 candidates, were created. From each set, the one pattern (out of 100) with the lowest Point Spread Function was selected. Hence, 300 patterns were produced for each undersampling ratio, making a total of 300 x 7 = 2100 patterns. The 2100 patterns were applied to a library of six fully sampled T 1 weighted images of Aplysia californica buccal ganglia. In order to compare the CS and fully sampled images, the l2-norm relative errors (RE) were calculated according to the equation:

√∑ ( ) √∑ (1)
where and are the signal intensities of voxel i in the fully encoded image and the undersampled image, respectively, and n is the total number of pixels. For each undersampling ratio, the averaged RE over the six images in the library was computed. The DLA undersampling pattern with the lowest RE was selected for each of the seven undersampling ratios and implemented in Paravision 5.1 (Bruker BioSpin, Ettlingen, Germany) starting from the standard FLASH pulse sequence.

Examples of k-space undersampling patterns are shown in Erreur ! Source du renvoi introuvable.. The k-space was undersampled along the two phase encoding directions and the pattern was repeated for every point in read direction.

Data acquisition

All MRI acquisitions were performed at 19°C on a 17.2 T system (Bruker BioSpin, Ettlingen, Germany) equipped with 1 T/m gradients. The RF transceivers used for imaging were custombuilt solenoidal microcoils with an inner diameter of 2.4 mm, the design of which has been described previously [START_REF] Jelescu | Highlighting manganese dynamics in the nervous system of Aplysia californica using MEMRI at ultra-high field[END_REF]. Two types of acquisitions were performed for each sample: a RARE respectively. Since the FOV size was not found to influence the DLA performance, the two groups were pooled together.

Sample preparation

A total of fourteen Aplysia californica were used in this study. Images acquired on ganglia from six animals were used for generating the library necessary to optimize the DLA based CS trajectories. Six other animals were used for acquiring fully sampled data sets. Besides providing reference images, these data sets were retrospectively undersampled in order to determine the optimal undersampling ratio. Finally, two animals were used to acquire both fully encoded and CS images. These data were also retrospectively undersampled. For all experiments the animals were food deprived for 48h prior to the beginning of the experiment in order to increase their food seeking behavior and maximize the intracellular Mn 2+ accumulation as described previously [START_REF] Radecki | Functional magnetic resonance microscopy at single-cell resolution in Aplysia californica[END_REF]. On the day of the experiment, the animals were injected with 100 mM MnCl 2 solution (500µl neurons, some of which are up to 200 µm in diameter [START_REF] Gardner | Bilateral Symmetry and Interneuronal Organization in the Buccal Ganglia of Aplysia[END_REF][START_REF] Gardner | Diphasic postsynaptic potential: a chemical synapse capable of mediating conjoint excitation and inhibition[END_REF] and can therefore be resolved with the spatial resolution employed here.

per
manner: the signal intensity of each voxel in the transverse plane (perpendicular to the longitudinal axis of the receiver coil) was normalized against the mean signal intensity of all voxels corresponding to ASW in this plane [START_REF] Radecki | Functional magnetic resonance microscopy at single-cell resolution in Aplysia californica[END_REF].

Results

Examples of fully encoded and retrospectively DLA-CS undersampled T 1 weighted images of the buccal ganglia are shown in Figure 3Erreur ! Source du renvoi introuvable..

In a first step we estimated the performance of the DLA-CS approach for seven different undersampling ratios both in terms of image resolution (characterized by PCC) and relative signal intensity error (Figure 4). As can be seen in Figure 4a, the PCC between the CS and fully encoded images for undersampling ratios higher than 50% drops to values below 0.8, generally considered as the threshold for a strong correlation [START_REF] Mukaka | Statistics Corner: A guide to appropriate use of Correlation coefficient in medical research[END_REF]. We observe an increase in PCC when averaging the signal over two repetitions. The relative signal intensity analysis was performed by comparing non-normalized fully encoded images with CS images retrospectively undersampled from the same raw data set. Measurements were performed in water (water ROI) and neuron bodies (cells ROI). For the latter the signal intensities were measured in the five biggest cells (B1, B2, B3, B6 and B9) and averaged. The relative error between the average signal intensities of fully encoded and CS images, for all the undersampling ratios considered, are displayed in Figure 4b. We notice that the error is inferior to 6% for all undersampling ratios considered.

However, for undersampling ratios larger than 60% the error corresponding to the cells ROI and the water ROIs diverge, which could introduce a bias in the signal intensity quantification. The different behavior of the signal intensity error in water versus cell bodies at large undersamplings is most likely due to the loss in spatial resolution (increased blurring) as indicated by the PCC results. Surprisingly, the signal intensity error did not show SNR dependence. However, not only does the standard deviation of the error for the cell ROIs increase with the undersampling ratio (Fig. 4b) but we also found it to be significantly higher for one repetition when compared to two repetitions (Supporting Table S1). Specifically, a student t-test showed a statistical significance (p = 0.0005) while for the water ROI no difference was found (p = 0.26).

Based on the results presented above we chose an undersampling of 50% for our next experiments which aimed at evaluating the performance of DLA-CS acquisitions for single neuron signal intensity quantification. Figure 5 presents Bland-Altman plots showing the difference in the normalized signal intensity values measured in single neurons in the fully sampled dataset and the 50% prospectively undersampled dataset (a) or the 50% retrospectively undersampled dataset (b). The average difference in the normalized signal intensity values was 1.3% and 0.5% for the prospectively and retrospectively undersampled data sets, respectively.

Conclusion and Discussion

It has been shown previously that MEMRI can be used to perform functional imaging of the buccal network of Aplysia californica as the manganese ions accumulate differentially in animals exposed to different food stimuli [START_REF] Radecki | Functional magnetic resonance microscopy at single-cell resolution in Aplysia californica[END_REF]. The purpose of this work was to reduce the acquisition time of such MEMRI protocols by implementing DLA-CS undersampling patterns ( 6) and evaluate whether a quantitative analysis of the signal enhancement remains feasible. To address this issue,

we obtained and compared fully encoded and DLA compressed T 1 weighted FLASH images of Aplysia californica buccal ganglia.

The performance of the DLA-CS FLASH acquisition was evaluated at various undersampling ratios. We found an undersampling ratio of 50% acceptable both in terms of image resolution and signal intensity quantification. Regarding single neuron signal intensity quantification we found, on average, a 1.37% percentage error between the fully sampled and prospectively undersampled data. This error was observed to be higher than the error measured using retrospective undersampling of the fully encoded data, which was found to be 0.50%. The difference between the two undersampling scenarios can be due to experimental errors such as hardware instability, subtle changes in the sample position in the B 0 field (resulting from vibrations associated with the strong encoding gradients), or slight sample modification.

When evaluating the Pearson Correlation Coefficient between the fully encoded and the undersampled images we notice that the performance of the DLA-CS technique is influenced by the image signal to noise ratio, in agreement with earlier studies (5), suggesting that higher accelerations are possible for higher SNR data. To summarize, our results demonstrate that the DLA-CS strategy proposed here can significantly accelerate data collection in high resolution quantitative MEMRI studies of neuronal tissues. Even though the acquisitions times remain long when compared to fast techniques, such as EPI or spiral imaging, the DLA-CS appears to be a promising approach at high magnetic fields and high spatial resolutions, where single shot acquisitions are not feasible. Moreover, the DLA-CS is not limited to magnetic resonance microscopy and could be also applied to preclinical and clinical studies, where shortening the acquisition time is equally desirable. Supporting Table S1: Signal intensity errors (Error) and corresponding standard deviations (SD) between fully encoded and compressed sensing acquisitions for one and two repetitions calculated in cells and water ROIs.

  acquisition, providing T 2 contrast (TR = 3000 ms, TE = 20 ms, A F = 4, 25 µm isotropic resolution) and a FLASH acquisition providing T 1 contrast (TR = 150 ms, TE = 2.441 ms, 3 averages, 2 repetitions, 25 µm isotropic resolution) in fully encoded and CS variants. The FOV was either 10x2.2x2.2 mm 3 or 10x2.0x2.0 mm 3 corresponding to matrix sizes of 400x88x88 and 400x80x80 and fully encoded FALSH acquisition times, per repetition, of 58 and 48 minutes,

  100g body weight; NaCl 345 mM, KCl 10mM, MgCl 2 25mM, MnCl 2 100mM, pH = 7.5) and were left in the aquarium for 45 minutes with unrestricted access to food (seaweed). The animals were then anesthetized with isotonic MgCl 2 solution (50ml per 100g body weight; MgCl 2 360mM, HEPES 10mM, pH = 7.5). Buccal ganglia were resected and inserted in 1.5 mm ID borosilicate glass capillaries (Vitro-Com, Mountain Lakes, NJ, USA) containing artificial sea water (ASW; NaCl 450 mM, KCl 10mM, MgCl 2 30mM, MgSO 4 20mM, pH = 7.5) and then slid inside the transceiver for MRI. The Aplysia buccal ganglia consists of approximately 300

Figure 1 .

 1 Figure 1. DLA acquisition patterns for 50%, 70% and 90% undersampling ratios. The horizontal axis represents the read direction, and the other two axes represent the phase encoding directions, as shown in the lower right corner.

Figure 2 .

 2 Figure 2. Schematic representation of the ROI selection for signal intensity quantification. Two acquisitions were performed for each sample: a T 2 weighted RARE (a), providing information about the sample anatomy and a T 1 weighted FLASH (b) reflecting the intracellularly accumulated Mn 2+ ions. Neurons were manually segmented on the RARE image and the corresponding ROIs were co-registered with the FLASH image. The drawn ROIs correspond to neurons B9 (red), B6 (orange) and B3 (green). Spatial resolution: 25 m isotropic.

Figure 3 .

 3 Figure 3. Fully encoded FLASH image (leftmost) and corresponding undersampled images after CS reconstruction (50%, 70% and 90% undersampling ratios). The CS images shown here were obtained by retrospectively undersampling the fully encoded k-space data.

Figure 4 .

 4 Figure 4. Evaluation of the performance of CS-DLA for different undersampling ratios. (a). PCC between fully and CS encoded images as a function of undersampling ratio for one (red) and two (blue) repetitions. The PCC was calculated according to Eq. 2 over the ganglia region containing approximately 50 000 voxels. Error bars represent standard deviations. (b) Percentage signal intensity difference between the fully encoded and CS images (the data represented is the average over two repetitions). Blue and red marks correspond to cell bodies and water regions, respectively. Error bars represent standard deviations. The data was obtained from 6 samples.

Figure 5 .

 5 Figure 5. Bland-Altman plots showing the difference in the normalized signal intensity values estimated from the fully sampled dataset and the 50% prospectively undersampled dataset (a) and the 50% retrospectively undersampled dataset (b). Each point corresponds to the signal intensity measured in one single neuron (2 samples, 10 neurons per sample).

Figure 1 Figure 2 Figure 3 13 Figure 4 Figure 5

 1231345 Figure 1

Acknowledgements

This work was funded by grant ANR-13-BSV5-0014-01 (project ANImE) and by the IRTELIS PhD program. Authors would like to acknowledge Boucif Djemai for assistance with animal husbandry.

Image reconstruction and analysis

The fully encoded images were processed directly in Paravision 5.1 software. CS undersampled data were reconstructed following the Split-Bregman algorithm provided by Goldstein and Osher in 2009 [START_REF] Goldstein | The Split Bregman Method for L1-Regularized Problems[END_REF], which was extended for this study to a 3D version with total variation (TV) penalty and Haar wavelet transform. Briefly, if we denote the undersampled Fourier transform corresponding to the undersampled pattern and is the Haar wavelet transform, the reconstructed image m is obtained by solving the following optimization problem:

Where y is the undersampled k-space data, µ is the expected noise level and

, where represents the value of gradient in k=x,y,z direction of image m at voxel index i. The original 2D code can be found in reference (13).

The performance of the DLA-CS strategy was evaluated by comparing signal intensities measured in individual neurons and in water in both fully encoded and undersampled images.

Five biggest neurons in the Aplysia's buccal ganglia (B1, B2, B3, B6 and B9) were identified and manually segmented on RARE (T 2 weighted) images (Figure 2a). (Note that as the buccal ganglia are bilaterally symmetric, one sample contains two neurons of each type.) The corresponding ROIs were co-registered to the FLASH (T 1 weighted) images and the mean signal intensity for each of them was calculated (Figure 2b).

In order to evaluate the extent of resolution loss between the fully encoded and undersampled images, we computed the Pearson's Correlation Coefficient (PCC) [START_REF] Bulusu | Determination of secondary flow morphologies by wavelet analysis in a curved artery model with physiological inflow[END_REF]:

Where and represent signal intensities of voxel i in the fully encoded and undersampled image, respectively, and and represent the mean signal intensity values of all voxels in the fully encoded and the undersampled image, respectively. ThePCC was calculated in manually drawn ROIs containing the ganglia and encompassing approximately 50 000 voxels.

The signal intensity quantification was performed by normalization against the water signal. To correct for possible RF inhomogeneities, the images were normalized in a position-dependent


Gabrielle Fournet

Jing-Rebecca Li

Denis Le Bihan

Dr Luisa Ciobanu
email: luisa.ciobanu@cea.fr

The influence of acquisition parameters on the metrics of the bi-exponential IVIM model

Keywords: IVIM, microvasculature, inflow effects, pulse sequence, diffusion encoding time IVIM, microvasculature, inflow effects, pulse sequence, diffusion encoding time

Purpose:

The IntraVoxel Incoherent Motion (IVIM) MRI signal, typically described as a monoexponential decay, can sometimes be better modeled as a bi-exponential function accounting for two vascular pools, capillaries and medium-size vessels. The goal of this work is to define precisely in which conditions the IVIM signal shape becomes bi-exponential and to understand the evolution of the IVIM outputs with different acquisition parameters.

Methods:

Rats were scanned at 7T and 11.7T using diffusion-weighted pulsed-gradient spinecho (SE) and stimulated-echo (STE) sequences with different repetition times (TR) and diffusion encoding times. The obtained IVIM signals were fit to the mono-and bi-exponential models and the output parameters compared.

Results:

The bi-exponential and mono-exponential models converge at long diffusion encoding times and long TRs. The STE is less sensitive to inflow effects present at short TRs, leading to a smaller volume fraction for the fast pool when compared to the SE sequence.

Conclusion:

The two vascular components are more easily separated at short diffusion encoding times, short TRs and when using SE sequences. The volume fractions of the two blood pools depend on the pulse sequence, TR and diffusion encoding time while the pseudodiffusion coefficients are only affected by the diffusion encoding time.

INTRODUCTION

Intravoxel Incoherent Motion (IVIM) magnetic resonance imaging (MRI) is getting momentum as an imaging modality capable of non-invasively extracting perfusion related parameters without the need for tracers or contrast agents. Initially, a mono-exponential model was proposed to describe the IVIM signal, taking into account only the smallest blood vessels, i.e. the capillaries [START_REF] Bihan | Separation of diffusion and perfusion in intravoxel incoherent motion MR imaging[END_REF]. To provide a more accurate description of the IVIM signal, we recently introduced (2) a biexponential IVIM model consisting of two components representing two distinct vascular pools.

One pool is associated with slow flowing blood in small vessels such as capillaries (as in the standard mono-exponential IVIM model) while the other pool corresponds to faster flowing blood in medium-size vessels, such as arterioles and venules (Erreur ! Source du renvoi introuvable.). In this manuscript, we establish under which experimental conditions the two pool model should be considered and study the evolution of its outputs with various acquisition parameters. In addition, we compare the IVIM outputs obtained with two different sequences often used in IVIM imaging: spin-echo and stimulated-echo. Ultimately, the results presented here will facilitate the interpretation of IVIM data acquired on different platforms and under different experimental protocols.

METHODS

Animals

Dark Agouti male rats (number of animals, n = 8, 232-349 g, 5-20 months, Janvier, Saint Isle,

MRI acquisitions

The acquisitions were performed on horizontally oriented 7T and 11.7T small animal MRI scanners (Biospec, Bruker Biospin, Etlingen, Germany) equipped with 760 mT/m gradient coil systems. 3 x 3 cm² four-element phased-array receiver coils and 7.2 cm (inside diameter) volume transmit coils (Bruker BioSpin, Etlingen, Germany) were used on both systems.

Coronal images were obtained using diffusion-weighted (Dw) pulsed-gradient spin-echo (PG-SE) or stimulated-echo (PG-STE) echo planar imaging (EPI) sequences. The common parameters to all trials were: 30 b-values (20 ranging from 5 to 500 s/mm² and 10 from 500 to 2500 s/mm², slice thickness 1.5 mm, field of view 20 x 20 mm 2 , 1 segment, 2 slices, diffusion gradient duration time  = 3 ms, gradient direction [X=0, Y=1, Z=0] (considering diffusion anisotropy is not present in the cortex at this resolution) and NR=6 repetitions. The other parameters were varied according to:

Data analysis

The images were processed using an in-house software written in MATLAB (MathWorks, Massachusetts, USA) as described in [START_REF] Fournet | A two pool model to describe the IVIM cerebral perfusion[END_REF]. Briefly, the data were averaged over the different repetitions in an ROI drawn manually on the cortical gray matter of the left hemisphere (shown in white in Erreur ! Source du renvoi introuvable..A).

The signal attenuation, ( ), including both the diffusion, ( ), and the IVIM, ( ), components, can be expressed as:

( ) ( ) ( ), [ 1
where is the diffusion coefficient of water in blood and the pseudo-diffusion coefficient and to the recently proposed bi-exponential IVIM model ( 2),

( ) ( ), [4]
where and , and , are the relative blood volume fractions and pseudodiffusion coefficients of the slow and fast pools, respectively. was set to 1.75 10 -3 mm²/s (4).

For experiments performed at 7T b lim was set to 500 s/mm 2 as in reference ( 2), however for the data acquired at 11.7T we used b lim = 300 s/mm 2 as we found that there was no IVIM signal left past this value probably due to the shortening of the T 2 of blood with the increase in the strength of the magnetic field.

Statistical methods

To determine the best model to fit the IVIM signal, the corrected Akaike information criterion, AICc, was calculated using MATLAB (

, [START_REF]Information theory and an extension of the maximum likelihood principle[END_REF] where N b is the number of b-values used to fit the signals, MSE the mean squared error and k the number of parameters in the model. Taking into account that the Gaussian noise hypothesis for the signal residuals counts as 1 parameter according to the AIC theory we obtain k = 3 and 5

for the mono-and bi-exponential models, respectively.

The difference in AICc between the two models, AICc mono -AICc bi , was used as estimation of the bi-versus the mono-exponential behavior of the IVIM signal; a positive difference indicating that the bi-exponential model performs better.

Using the R software (R Foundation for Statistical Computing, Vienna, Austria) (6), Welch's ttests [START_REF] Welch | The generalization of "Student's" problem when several different population variances are involved[END_REF] were performed to assess whether the output parameter means were significantly different from each other while varying the acquisition parameters. A P-value < 0.05 was considered statistically significant. The minimum sample size to perform the statistical tests was also estimated and established to be respected for all IVIM parameters found significant in this study. At TRs of 3500 ms, the difference in AICc becomes negative and the IVIM signal is better fit to a mono-exponential model for both sequences. There is no significant difference in f IVIM , D* and the difference in AICc between the two pulse sequences at this value of TR.

RESULTS

(c) Influence of the diffusion encoding time

The STE sequence was employed with three different values: 14, 30 and 60 ms. The IVIM signal is better fit to the bi-exponential IVIM model for all three -values, but f IVIM and the difference in AICc between the two IVIM models are not significantly different between the different -values although the latter decreases with increasing . f fast significantly decreases between  = 30 and 60 ms. D* slow is significantly smaller for = 14 ms compared to both  = 30 and 60 ms (Table 3). Finally, D* fast is shown to significantly decrease with increasing -value.

DISCUSSION

By decreasing TR, we observe an increase in f IVIM . This increase is coherent with inflow effects which are present at short TRs in 2D sequences employing a small number of slices [START_REF] Gao | Inflow effects on functional MRI[END_REF]. The volume fraction of the slow pool, f IVIM *f slow , does not vary significantly with TR while the volume fraction of the fast pool, f IVIM *f fast , is significantly larger at short TRs. This is consistent with inflow effects being more important for faster flows. Note that the inflow effects will be significantly reduced in studies employing multi-slice or 3D acquisitions.

The bi-exponential behavior of the IVIM signal seems to be less visible at long TRs as the difference in AICc between the two IVIM models decreases. This implies that, at longer TRs, the two models converge, most likely because the fast flowing pool is harder to detect than at short TRs where its contribution is enhanced by inflow effects.

At short TRs, the two pulse sequences used, SE and STE, present different sensitivities to inflow effects with the STE sequence being less affected as the values measured for f IVIM and f fast are significantly lower than for the SE sequence. This is due to the fact that the signal from inflowing spins is decreased in the STE compared to the SE sequence due to the higher first moments of the slice-selection gradients [START_REF] Caprihan | Velocity imaging of slow coherent flows using stimulated echoes[END_REF]. This differential sensitivity to inflow is further confirmed by the fact that f IVIM *f slow does not differ significantly with the pulse sequence whereas f IVIM *f fast decreases significantly between the SE and STE sequences. Additionally, as the fast pool, exhibiting stronger inflow effects, is less visible with the STE sequence the difference in AICc between the two IVIM models is significantly lower for the STE sequence than for the SE sequence.

At long TRs for which inflow effects are negligible, the IVIM signal becomes mono-exponential for the two pulse sequences. This was predictable from the evolution of the IVIM signal with TR as discussed in the previous paragraph. Without inflow effects, the IVIM parameters are not significantly different for the two pulse sequences.

The dependence of the IVIM parameters on the diffusion encoding time has already been partly discussed previously [START_REF] Fournet | A two pool model to describe the IVIM cerebral perfusion[END_REF]. To further be able to increase the diffusion encoding time without increasing the echo time, the STE sequence was used instead of the SE sequence in this study.

As observed previously, a decrease in the difference in AICc is observed when increasing  suggesting that the two IVIM models, mono-and bi-exponential, converge at long diffusion encoding times. Accordingly, the fraction of fast flowing blood vessels decreases for the longest . The pseudo-diffusion coefficients of the two pools are also affected. D* slow increases between  = 14 and 30 ms implying that the slow pool is likely to be in the sinc regime at the shortest  for which D* can be expressed as D* = V²/6 (where V is the mean blood velocity). As , increases, D* slow reaches the exponential regime for which D* = LV/6 (where L is the mean vessel length). In this regime, D* is independent of . On the contrary, D* fast was found to decrease with . Indeed, as  increases, the fastest spins cannot be detected anymore leading

CONCLUSIONS

This study confirms the contribution to the IVIM signal of flow within both small vessels (capillaries) and larger vessels (on the scale of arterioles and venules). We have shown that, in the case of 2D spin-echo acquisitions with a small number of slices, the main factor influencing the behavior, bi-or mono-exponential, of the IVIM signal as well as the measured volume fractions is the repetition time. Shorter repetition times allow for an easier separation of the two pools. However, they can lead to an overestimation of f IVIM and f fast. As expected, the IVIM output parameters depend less on the repetition time for sequences less sensitive to inflow effects such as the stimulated echo sequence. The diffusion encoding time also influences the IVIM metrics.

For very long diffusion encoding times, the behavior tends to become mono-exponential. For certain ranges, this parameter can also influence the values obtained for the pseudo-diffusion coefficients of the slow and fast pools.

To conclude, certain acquisition parameters can have an effect on the IVIM outputs and therefore care should be taken when comparing results obtained under different experimental protocols. Tables:

Table 1. IVIM parameters for the SE sequence and two TRs (mean ± SD, n = 4

  France) were used in this study. All animal experiments were conducted according to recommendations of the EU Directive 2010/63/EU for care and use of laboratory animals. The protocol was approved by the Comité d'EThique en Expérimentation Animale Commissariat à l'Energie Atomique et aux énergies alternatives Direction des Sciences du Vivant Ile de France (CETEA CEA DSV IdF) under protocol IDs 10_032 and 15_040.

  (a) For investigating the influence of TR (experiments performed at 7T): Dw-PG-SE-EPI sequence with TE = 45 ms, TR = 1000 and 3000 ms, diffusion gradient separation time  = 14 ms, in-plane resolution 250 x 250 µm², 6 averages and n=4. (b) For investigating the influence of the pulse sequence (experiments performed at 7T): Dw-PG-SE-EPI and Dw-PG-STE-EPI sequences with TE= 45ms , TR=1000 and 3500 ms,  = 14 ms, in-plane resolution 400 x 400 µm², matrix size 50 x 50, 6 averages and n=6 and 4 for the two TRs, respectively. (c) For investigating the influence of the diffusion encoding time (experiments performed at 11.7T): Dw-PG-STE-EPI sequence with TE/TR = 18/1000,  = 14, 30 and 60 ms, inplane resolution 400 x 400 µm², 4 averages and n = 6.

  and the overall signal when b=0. The diffusion parameters were first obtained by fitting ( ) for b-values > b lim to the diffusion Kurtosis model (3), diffusion coefficient obtained when b approaches 0 and is the Kurtosis parameter characterizing the deviation from the exponential decay. The diffusion component was extrapolated for b-values < b lim and subtracted from ( ). The remaining IVIM component was fit to the standard mono-exponential IVIM model (1),

( a )

 a Influence of TRErreur ! Source du renvoi introuvable. gathers the results obtained for TR = 1000 ms and 3000 ms. Both f IVIM and f fast show a significant increase with decreasing TR, with f IVIM *f fast significantly higher at short TRs. The other IVIM parameters do not vary significantly. The difference in AICc decreases with increasing TR, suggesting that the bi-exponential behavior is more present at short TRs.(b) Influence of the pulse sequenceThe two pulse sequences, SE and STE, were compared for two TR values, 1000 and 3500 ms.The results are displayed in Erreur ! Source du renvoi introuvable.. At short TRs, f IVIM , f fast and f IVIM *f fast are significantly smaller for the STE sequence. Although not significantly different between the two sequences, the difference in AICc is smaller for STE, suggesting that the biexponential behavior is less visible with this sequence compared to the SE sequence.

  to a decrease in D* fast .Erreur ! Source du renvoi introuvable. summarizes the impact of different combinations of acquisition parameters studied in this work on the IVIM signal behavior and the IVIM outputs. At short TRs, with both pulse sequences, the IVIM signal behavior is bi-exponential with f IVIM and f fast overestimated for the SE sequence at all diffusion encoding times and for the STE sequence only at short diffusion encoding times. At long TRs, the mono-and bi-exponential models converge, with the mono-exponential model being a better fit at long diffusion encoding times. D* slow and D* fast are only impacted by the diffusion encoding time.When inflow effects are present the signal contribution from blood in the fast pool is artificially increased. One should also consider that the flowing blood fractions of the IVIM model (f IVIM , f fast and f slow ) are T 1 and T 2 weighted. Hence, in order to estimate the correct values for f fast , f slow and f IVIM from the measured values one should take into account the relaxation times, T 1 and T 2, of blood and tissue. While the influence of the echo time was not investigated in this study, it is expected to be significant for the estimation of the volume fractions of the two pools. An extended bi-exponential IVIM model correcting for the TE and TR dependences, paralleling the extended mono-exponential model recently proposed and applied to other organs[START_REF] Jerome | Extended T2-IVIM model for correction of TE dependence of pseudo-diffusion volume fraction in clinical diffusion-weighted magnetic resonance imaging[END_REF], will certainly increase the reliability of the IVIM measurements.

Figure 1 .

 1 Figure 1. The bi-exponential model reflects two vascular pools: capillaries and larger vessels. (A) MRI image showing a coronal section of the rat brain. (B) Drawing of the architecture of the microvasculature consisting of arterioles, venules and capillaries.

Table 2 .

 2 ). P-values < 0.05 are highlighted in bold. IVIM parameters for SE and STE pulse sequences at two TRs (mean ± SD). P-values < 0.05 are highlighted in bold.

		Output parameters TR = 1000 ms TR = 3000 ms		P-value
		f IVIM (%)	13.41 ± 0.55		6.36 ± 1.06	0.000076
		f fast (%)	75.61 ± 10.84 42.06 ± 10.36		0.0021
		f IVIM *f slow (%)	3.25 ± 1.45		3.76 ± 0.92		0.29
		f IVIM *f fast (%) D* slow (10 -3 mm²/s) D* fast (10 -3 mm²/s)	10.16 ± 1.51 2.21 ± 0.33 27.48 ± 1.97		2.60 ± 0.77 2.16 ± 0.70 25.96 ± 3.21	0.00026 0.45 0.23
		AICc mono -AICc bi	19.04 ± 13.18	4.35 ± 7.57		0.057
		TR = 1000 ms (n = 6)			TR = 3500 ms (n = 4)
	Output Parameters	SE	STE		P-value	SE		STE	P-value
	f IVIM (%)	10.84 ± 1.16	6.34 ± 1.88	0.00046	3.41 ± 0.56	3.27 ± 1.19	0.42
	f fast (%)	59.50 ± 14.29	29.98 ± 19.88		0.0080	-		-	-
	f IVIM *f slow (%)	4.42 ± 1.63	4.67± 1.92		0.41	-		-	-
	f IVIM *f fast (%)	6.42 ± 1.69	1.67 ± 1.21	0.00016	-		-	-
	D* slow (10 -3 mm²/s)	2.37 ± 1.62	2.37 ± 0.13		0.50	-		-	-
	D* fast (10 -3 mm²/s)	21.45 ± 4.19	17.02 ± 5.04		0.065	-		-	-
	D* (10 -3 mm²/s)	-	-		-	6.37 ± 1.19	6.72 ± 2.64	0.41

ACKNOWLEDGMENTS

The authors thank Brad Sutton and Alex Cerjanic for useful discussions and Boucif Djemai and Erwan Selingue for their excellent technical assistance in the animal experiments. The research was supported by the ANR Grant ANR-13-NEUC-0002-01.

Note. The symbols =, -, + and ++ signify, respectively, equal, smaller, higher or much higher values compared to those obtained for long TR and short .