Olivia Miraucourt

Stéphanie Salmon
email: stephanie.salmon@univ-reims.fr

Marcela Szopos

Marc Thiriet

Blood flow in the cerebral venous system: modeling and simulation

Keywords: intracranial venous flow, computational fluid dynamics, Navier-Stokes equations, 3D reconstruction

The development of a software platform incorporating all aspects, from medical imaging data, through three-dimensional reconstruction and suitable meshing, up to simulation of blood flow in patient-specific geometries, is a crucial challenge in biomedical engineering. In the present study, a fully three-dimensional blood flow simulation is carried out through a complete rigid macrovascular circuit, namely the intracranial venous network, instead of a reduced order simulation and partial vascular network. The biomechanical modeling step is carefully analyzed and leads to the description of the flow governed by the dimensionless Navier-Stokes equations for an incompressible viscous fluid. The equations are then numerically solved with a free finite element software using five meshes of a realistic geometry obtained from medical images to prove the feasibility of the pipeline. Some features of the intracranial venous circuit in the supine position such as asymmetric behavior in merging regions are discussed.

Introduction

The cerebral vasculature is a highly complicated three-dimensional network, the role of which is the maintenance of a proper blood drainage. The large cerebral arteries (e.g. carotid arteries, basilar trunk, and the circle of Willis) are the focus of multiple investigations [START_REF] Formaggia | Basic mathematical models and motivations[END_REF] and references therein). In this context, blood flow computational models are relevant to various clinical applications, such as check-up and follow-up of aneurysms [START_REF] Cebral | Efficient pipeline for imagebased patient-specific analysis of cerebral aneurysm hemodynamics: Technique and sensitivity[END_REF]; [START_REF] Shojima | Computational fluid dynamic simulation of a giant basilar tip aneurysm with eventual rupture after hunterian ligation[END_REF]; [START_REF] Retarekar | Stratification of a population of intracranial aneurysms using blood flow metrics[END_REF]) and control after stenting [START_REF] Wang | Computational haemodynamics in two idealised cerebral wide-necked aneurysms after stent placement[END_REF] or coiling [START_REF] Morales | A virtual coiling technique for image-based aneurysm models by dynamic path planning[END_REF]). In contrast, the behavior of the venous cerebral system is rarely explored. This compartment has its own features, such as the lack of parallelism between arterial and venous circulations [START_REF] Schaller | Physiology of cerebral venous blood flow: From experimental data in animals to normal function in humans[END_REF]) and asymmetrical and much more varying between-subject architecture and behavior compared to the arterial circuit [START_REF] Schaller | Physiology of cerebral venous blood flow: From experimental data in animals to normal function in humans[END_REF]; [START_REF] Stoquart-Elsankari | A phasecontrast MRI study of physiologic cerebral venous flow[END_REF]).

In parallel to technological improvements of medical imaging techniques [START_REF] Ozsarlak | MR angiography of the intracranial vessels: Technical aspects and clinical applications[END_REF]; [START_REF] Miyazaki | Nonenhanced MR angiography[END_REF]), computational fluid dynamics evolved with the development of new numerical schemes and high performance computing. However, despite recent progress in blood flow simulations, open problems still exist. They derive not only from computational challenges and cost issues, but also from the difficulty to accurately model the multiple physiological assumptions of these complex multiphysics and multiscale phenomena. In this methodological framework, two strategies are generally considered.

The first strategy consists in carrying out reduced order simulations, i.e., to consider either zero-dimensional (0D) lumped-parameter models, exploiting the electric circuit analogy for the circulatory system, or one-dimensional (1D) distributed parameter models, thereby dealing with large vascular networks: [START_REF] Ursino | A simple mathematical model of the interaction between intracranial pressure and cerebral hemodynamics[END_REF] (0D model of intracranial hemodynamics), [START_REF] Gadda | A new hemodynamic model for the study of cerebral venous outflow[END_REF] (lumped parameter model of the cerebral venous outflow), [START_REF] Reymond | Validation of a one-dimensional model of the systemic arterial tree[END_REF] (1D model of the arterial network), [START_REF] Ho | Numerical simulation of blood flow in an anatomically-accurate cerebral venous tree[END_REF] (1D flow model built on an anatomical geometry of the cerebral venous system) and Müller, Toro (2014a), Müller, Toro (2014b) (closed-loop 1D model of the human circulation, with a special focus on the venous system).

Alternatively, three-dimensional (3D) computational fluid dynamics simulations can be carried out in patient-specific geometries. Classically, modeling is limited to relatively small samples of the vascular circuit, possibly coupled to 0D and 1D models, whether or not wall deformability is incorporated in modeling [START_REF] Formaggia | Basic mathematical models and motivations[END_REF]; [START_REF] Cebral | Efficient pipeline for imagebased patient-specific analysis of cerebral aneurysm hemodynamics: Technique and sensitivity[END_REF]; [START_REF] Boissonnat | From arteriographies to computational flow in saccular aneurisms: The INRIA experience[END_REF]). Recent attempts were proposed for the cerebral arterial network: [START_REF] Passerini | A 3D/1D geometrical multiscale model of cerebral vasculature[END_REF] (1D description of the circle of Willis, coupled to a fully 3D model of a carotid artery), [START_REF] Xiao | Multi-scale computational model of three-dimensional hemodynamics within a deformable full-body arterial network[END_REF] (3D/0D computational model of hemodynamics within a deformable full-body arterial network), and venous system: [START_REF] Caiazzo | Computational haemodynamics in stenotic internal jugular veins[END_REF] (3D description of blood flow in rigid jugular veins coupled with a 1D model for major intracranial veins).

Our contribution is to rely on 3D blood flow modeling in more complete vascular networks, in order to better describe blood flow patterns in complicated geometries. In particular, such approach is relevant in the cerebral venous network whose geometry significantly differs from arteries. The present work is aimed at proposing a complete computational modeling approach from medical image processing to numerical resolution. The study relies on developing and using open-source softwares, which enable a full handling of the underlying methods and afford guarantees in terms of availability as well as reproducibility. We focus on intracranial venous blood flow in supine position, that is, a drainage through the jugular veins, which are nevertheless not included in the present model, as they are strongly deformable. In particular, the present model assumes that intracranial veins located between the skull and the brain are rigid. The venous circuit is isolated, as reduced modeling approaches (e.g., 0D and 1D models) that can be coupled to 3D models need patient-specific geometrical and mechanical data that cannot be easily assessed.

The remainder of this article is structured as follows. Section 2 is devoted to the description of the venous network, Section 3 to the biomechanical model and its associated equation set, and Section 4 to the numerical procedure. Results and discussion are presented in Sections 5 and 6. Concluding remarks and future prospects are gathered in the final section.

Anatomical model of the large cerebral veins

The large cerebral veins are localized in the leptomeninges and drain blood to the internal jugular veins that ensure the extracranial drainage of the brain perfusion. The realistic three-dimensional description of the geometry of the network was constructed starting from 3D angiographic images obtained by magnetic resonance angiography (MRA). More details on the different steps of the process can be found in [START_REF] Miraucourt | 3D CFD in complex vascular systems: A case study[END_REF].

The brain venous network is composed by input veins (Fig. 1): vein of Galen (7); internal cerebral vein (8); basilar vein (9); superior cerebral veins (10); superior anastomotic veins (11); draining the blood into the superior sagittal sinus (2) and straight sinus (3), down to their confluence (4). The blood then passes into the transverse (5) and sigmoid parts (6) of the lateral sinuses, and reaches their corresponding extracranial regions.

Depending on the subject's position (i.e., lying vs. standing), the cerebral blood flows back to the cardiac pump using different set of cerebral veins [START_REF] Gisolf | Human cerebral venous outflow pathway depends on posture and central venous pressure[END_REF]; [START_REF] Stoquart-Elsankari | A phasecontrast MRI study of physiologic cerebral venous flow[END_REF]). In the supine position, the functional drainage from the left and right side of the venous circuit encompass large jugular veins that can be included in the 3D reconstruction. In the standing position, a collection of small veins drain the cerebral blood. However, current medical image acquisition processes do not include a more sophisticated description of the collateral pathways (the vertebral plexus and anastomoses leading blood to the downstream segments of the jugular veins) that were very recently shown to play an important role as an alternative route for venous outflow according to the subject position, that is, standing or sitting [START_REF] Gadda | A new hemodynamic model for the study of cerebral venous outflow[END_REF]).

Neither the jugular veins, which drain the intracranial venous blood in the supine position, nor the collection of small veins, which drain the cerebral blood in the standing position, are included in the present model. However, two fictitious exit veins (Fig. 1 (1)) have been added to the computational domain, so that the outlet boundary condition does not affect the flow dynamics within the computational domain. Moreover, the downstream influence of a 3D unsteady flow can be physically mimicked. It is worth mentioning that in the considered venous network, valves are absent [START_REF] Kiliç | Anatomy of cerebral veins and sinuses[END_REF]), at the difference of the venous lower limb system, where a complex system of valves ensures the passage of blood upward from the feet against gravity.

Approximate values of vein radii, computed from the geometric vessel modeling process are compared with literature data in Table 1, with an emphasis on values previously used in 1D models of venous blood flow. These results are in good accordance, considering the individual variations.

Biomechanical model

The following standard blood properties and flow features and assumptions are used: (1) blood is incompressible;

(2) blood is homogeneous, as the size of red blood cells (97% of the flowing cell population) is much lower than the vein caliber;

(3) the endothelial wetted surface is smooth at this length scale (the size of indentations linked to endothelial cell nuclei is of the same order of magnitude than of that of blood cells); (4) blood is Newtonian; (5) venous walls are rigid;

(6) flow is isothermal; and (7) body forces are neglected (gravity and electromagnetic field as flow is computed at a time immediately after scanning, before standing up).

Rheology of blood and vascular walls still remain important modeling issues and challenging numerical tasks, that are not yet fully handled. The effect of stress history and blood thixotropy is, in general, not incorporated in blood flow models of the literature, thereby not representing properly blood rheology. Wall mechanics is represented, in general, by a simple shell model. In the absence of blood stagnation in the large compartment of the intracranial venous circuit (Table 1), that is, in the absence of red blood cell rouleau formation, blood can be supposed to remain Newtonian [START_REF] Formaggia | Basic mathematical models and motivations[END_REF]; [START_REF] Thiriet | Cell and Tissue Organization in the Circulatory and Ventilatory Systems[END_REF]). Hence, thixotropic, viscoelastic, and shear-thinning features are omitted.

Intracranial veins run through the meninges, a quasi-incompressible biological tissue composed mostly of water, which are quite constrained between the quasi-incompressible brain and rigid skull. They are thus not susceptible to experience large deformations at normal transmural pressures. Moreover, recent computational sensitivity analyses [START_REF] Sforza | Hemodynamic analysis of intracranial aneurysms with moving parent arteries: Basilar tip aneurysms[END_REF]) showed that neglecting wall motion is an acceptable hypothesis when studying intracranial arterial blood flow, although arteries are subjected to stronger changes in internal pressure. Wall rigidity assumption is more questionable for the jugular veins that drain the cerebral blood in the supine position, as these veins travel through the neck and are thus exposed to the ambient pressure and deformation engendered by pressure waves running along the carotid arteries (Section 6). In order to ensure wall rigidity hypothesis, the jugular veins are not included in the present model.

October 21, 2016 a journal

MirSalSzoThi

Dimensional analysis

As a consequence of all these assumptions, we consider in the following the governing equations of an unsteady isothermal flow of a Newtonian incompressible fluid (mass density ρ, dynamic viscosity µ, and kinematic viscosity ν = µ/ρ) in an anatomical vessel conveyed with a velocity u(x, t) (x: Eulerian position, t: time), derived from the mass and momentum conservation. They are known as the non-linear Navier-Stokes equations

ρ(∂ t + u • ∇)u = f + ∇ • C (1) ∇ • u = 0
where C = -pI + 2µD is the stress tensor, p is the pressure, I the metric tensor and

D = 1 2 (∇u + (∇u) T
) is the rate of deformation tensor for a Newtonian fluid. In isothermal conditions, Equations (1) lead to the simplified form of the Navier-Stokes equations in the absence of body forces f , namely

ρ(∂ t u + (u • ∇)u) = -∇p + µ∆u.
The Navier-Stokes equations can be transformed into a dimensionless form, using appropriate scales for the length (L ; either the boundary layer thickness or the duct hydraulic radius for fullydeveloped flow), time (T ; either the flow cycle period [T ≡ T ] or the inverse of the flow pulsation [T ∝ ω -1 ]), velocity (U ; either the mean or maximal cross-sectional velocity), and the pressure (P ). Therefore, we have t = t/T , x = x/L , ũ = u/U , p = p/P , and the dimensionless differential operators are then expressed by ∇ = L ∇ and ∆ = L 2 ∆.

The classical formulation of the dimensionless Navier-Stokes equations are hence obtained

ρU T ∂ t ũ + ρU 2 L (ũ • ∇)ũ = - P L ∇p + µU L 2 ∆ũ.
According to the predominant force type, either convective inertial forces or viscous forces, the dimensionless coefficients can be expressed differently. In steady flows, the pressure gradient balances the friction (and gravity when it has significant effect, as in blood vessels in standing position). In quasi-steady flows, the rate of change in boundary conditions are so slow that the momentum has time to diffuse during the flow period. In unsteady flows, the local inertia effects can be preponderant.

Since we are interested in the case when convective inertia is predominant, the dimensionless Navier-Stokes equations become (when tilde notation is omitted for the sake of simplicity)

St ∂ t u + (u • ∇)u = -∇p + Re -1 ∆u where Re = ρL * V * µ , St = ωL * V * =

Initial and boundary conditions

The initial status of the fluid velocity, u 0 has to be carefully prescribed, since it has to be divergence-free to be admissible. Unfortunately, in hemodynamic computations, u 0 is usually unknown, hence chosen equal to zero everywhere or, as a better guess, as the solution of a stationary Stokes problem.

The inflow is situated at the microcirculation exit, far from the thorax, hence far from the effects of breathing and cardiac pumping [START_REF] Thiriet | Biology and Mechanics of Blood Flows, part I: Biology of Blood Flows, part II: Mechanics and Medical Aspects of Blood Flows[END_REF]). Since flow is not fully developed, both Poiseuille (steady case) or Womersley (unsteady case) flows are not appropriate. An impulsive start, with a uniform injection velocity is then employed. However, rigorous analysis of the unsteady Navier-Stokes equations demands regularity of the flow solution down to t = 0, so this condition must be modeled as smoothly as possible. In practice, a time-dependent multiplier of the form (1 -e -kt ), with k constant, is introduced to ensure that the inflow data smoothly ramps up from zero to the desired steady state boundary condition (a constant function u in ) as t → ∞. This velocity inlet condition is imposed at the 29 entrance sections of the mesh. The value of u in is chosen in accordance to Table 2 and the parameter k is chosen such that the steady state value is reached at the third of the first cycle.

At the level of the exit veins, in a first approximation, we use a traction-free boundary condition and focus on the intracranial hemodynamics to minimize the influence of the outflow section. For more realistic boundary conditions, we report to another work carried out in the VIVABRAIN project1 [START_REF] Chabannes | Hemodynamic simulations in the cerebral venous network: A study on the influence of different modeling assumptions[END_REF]).

In summary, the imposed boundary conditions are the following: (i) a smoothed impulsive velocity profile at entries; (ii) homogeneous natural conditions at exits; and (iii) no-slip condition at walls assumed to be rigid.

Final model

The mathematical model that we consider in the sequel is the Navier-Stokes equations with mixed boundary conditions:

             ∂ t u + u • ∇u -ν∆u + ∇p = f in Ω × (0, T ) ∇ • u = 0 in Ω × (0, T ) u = g 1 on Γ D × (0, T ) ν ∂u ∂n -pn = g 2 on Γ N × (0, T ) u| t=0 = u 0 in Ω (2)
where Ω ⊂ R 3 is a bounded domain with a Lipschitz-continuous boundary Γ, (0, T ) is a finite time interval. We recall that u = (u 1 , u 2 , u 3 ) is the fluid velocity, p its pressure and ν its kinematic viscosity. We suppose that Γ consists of two measurable parts: Γ D , where Dirichlet boundary conditions are imposed, and Γ N , where Neumann boundary conditions are prescribed. We denote by n the unit outward normal vector to Ω on its boundary Γ. The Dirichlet boundary conditions correspond in the non-homogeneous case to a velocity g 1 , generally prescribed at the inlets, and in the homogeneous case to the non-slip boundary condition prescribed at the vessel wall; the function g 2 corresponds to the normal component of the stress tensor, prescribed at the outlets, which in the current model is set to zero.

Numerical method

As derived in the previous section, the equations under consideration are the unsteady incompressible Navier-Stokes equations (2), endowed with appropriate initial and boundary conditions. In this section, we describe the numerical approach we use in the sequel to compute an approximate solution to these equations, involving a time-scheme based on the characteristics method and a spatial discretization of finite element type. This algorithm is quite classical and is used as a reference to prove feasibility of the pipeline.

The main idea of the time discretization scheme is to take into account the fact that the terms ∂ t u + u.∇u can be written as the total derivative of u in the direction of the flow u [START_REF] Pironneau | On the transport-diffusion algorithm and its applications to the Navier-Stokes equations[END_REF]; [START_REF] Formaggia | Basic mathematical models and motivations[END_REF]). More precisely, let x ∈ Ω and t ∈ [0, T ]; we start by defining the characteristics X = X (τ ; t, x) associated to the vector field u as the solution to

d dτ [X (τ ; t, x)] = u(X (τ ; t, x), t), τ ∈ (0, T ) X (t; t, x) = x (3)
The total derivative can then be written as:

(∂ t u + u.∇u) (x, t) = d dτ [u(X (τ ; t, x), t)] | τ =t (4)
where X (τ ; t, x) is the position at time τ of a particle which has the position x at time t. Consequently, it can be approximated by:

(∂ t u + u • ∇u) (x, t n+1 ) = u(x, t n+1 ) -u(X n (x), t n ) ∆t + O(∆t) (5)
where

X n (x) = x -u(x, t n )∆t + O(∆t 2
) is an approximation of X (t n+1 ; t n , x), ∆t > 0 is the time-step and t n = n∆t, n ∈ N. As a consequence, the non-linear term is treated explicitly. The term containing X n (x) goes to the right hand side, which yields a generalized Stokes problem, that can be solved as a standard Stokes problem, and can be handled numerically as such, with a preconditioned Uzawa/conjugate gradient algorithm.

In order to present the space discretization method, we introduce the functional spaces for the velocity V = {v ∈ H 1 (Ω) | v |Γ D = 0} and for the pressure Q = L2 (Ω), respectively. Let T be a triangulation of the domain Ω and h the maximum of the diameters of the triangles of T . We introduce the finite-dimensional spaces V h and Q h which are respectively contained in V and Q 6)

V h = {v h ∈ C 0 ( Ω) d , ∀K ∈ T , v h|K ∈ P 2 (K)} (
Q h = {q h ∈ C 0 ( Ω), ∀K ∈ T , q h|K ∈ P 1 (K)} ( 7
)
where P 1 (K) (resp. P 2 (K)) is the space of polynomials of degree less than or equal to 1 (resp. 2) in each triangle of T . The discretization spaces are chosen compatible -in the sense of the so-called inf-sup condition -to cope with the divergence constraint (Formaggia et al. 2009, Chapter 15).

The fully discrete problem then reads: find

u h ∈ V h and p h ∈ Q h such that      1 ∆t u n+1 h , v h + ν ∇u n+1 h , ∇v h -p n+1 h , ∇ • v h = f n+1 h , v h + 1 ∆t (u n h • X n , v h ) + g n+1 2h , v h ΓN ∇ • u n+1 h , q h = 0 with u n+1 h = g n+1 1h on Γ D .
Theoretically expected convergence in space for velocity error is of order 3 for the L 2 -norm and of order 2 for the semi-norm H 1 ; for the pressure error, it is of order 2 for the L 2 -norm and of order 1 for the semi-norm H 1 [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations[END_REF], Chapter 2). The characteristics method is proven to be unconditionally stable and expected to be of order 1 in time (Formaggia et al. 2009, Chapter 15).

The numerical solving approach relies on the finite element library FreeFem++ 2 (Hecht (2012)) as described in (Miraucourt, Garnotel (2014)). To test the whole chain from images to computation, we developed an in-house software, based on the Visualisation ToolKit3 in the framework of the VIVABRAIN project. Several meshes are obtained from the 3D reconstruction (Table 4). We used the algorithm to compute flow with a time-step equal to 10 -2 s, thus providing a simulation of 2 seconds of a quasi-steady blood flow in the cerebral venous network. The CPU time with respect to the different meshes are summarized in field computed on a finite-element basis are visualized with color-coded magnitudes according to ParaView software4 and flow rates are calculated across selected cross-sections.

Results

We present in Fig. 2 the velocity field in the whole intracranial cerebral venous network. A complicated three-dimensional flow behavior is observed, especially in merging regions of the intracranial venous circuit. This asymmetric behavior appears for instance in a recirculation in the entrance segment of the right transverse sinus. The asymmetric feature of the venous flow can be explained, at least partly, by the asymmetric architecture of the venous network. The more refined the mesh, the clearer the recirculation.

Note that for comparison purposes, only MRI measurements in the supine position are extracted from the literature (which deals with geometrically reduced 1D models). In the present model, the predicted values of the velocity magnitude in the superior sagittal sinus range from 12.5 • 10 -2 m/s to 18 • 10 -2 m/s, consistently with clinical data, that is, 15 • 10 -2 m/s [START_REF] Schaller | Physiology of cerebral venous blood flow: From experimental data in animals to normal function in humans[END_REF]) and 15.2 ± 3 • 10 -2 m/s [START_REF] Gideon | Measurement of blood flow in the superior sagittal sinus in healthy volunteers, and in patients with normal pressure hydrocephalus and idiopathic intracranial hypertension with phase-contrast cine MR imaging[END_REF]) from a set of 14 control subjects using MR velocimetry. Flow rates calculated at inlets and outlets are identical (∼ 8 ml/s for an entry velocity u in = 50 mm/s, corresponding to a physiological value [START_REF] Stoquart-Elsankari | A phasecontrast MRI study of physiologic cerebral venous flow[END_REF]), hence guaranteeing mass conservation. Similar or slightly higher values, between 10 ml/s and 14 ml/s for the outflow were reported in the literature (Müller, Toro (2014a,b)). Comparisons between the model predictions of flow rate values at selected cross-sections with MRI measurements reported in the literature are presented in Fig. 3.

Discussion

Reduced dimensions models

Previous studies mainly focused on the arterial compartment or investigated the venous intracranial dynamics by means of reduced 0D or 1D models. However, local values of vascular resistances and compliances remain, in general, unknown in the context of personalized studies. Moreover, they cannot yield detailed flow patterns. In addition, it has been recently demonstrated in [START_REF] Xiao | A systematic comparison between 1-D and 3-D hemodynamics in compliant arterial models[END_REF] that a satisfactory agreement is obtained by 1D models during decelerating phase, but significant differences appear during accelerating phase, partly explained by the inability of current 1D models to take into account the helical motion induced by the vessel curvature. In the intracranial arterial network, a good agreement between both models can be obtained when arterial wall elasticity is ignored [START_REF] Grinberg | Modeling blood flow circulation in intracranial arterial networks: A comparative 3D/1D simulation study[END_REF]). However, fluctuations in the pressure and flow rates predicted by the rigid wall model are significantly higher than those predicted by the flexible wall model.

The present investigation provides information on blood flow behaviour using a three-dimensional personalized geometry. The 3D computational model showed that a jet can be observed during the beginning of the decelerating phase after peak flow in the entrance segment of the right transverse sinus downstream from a slight constriction, probably due to the printing of an adjoining anatomical element in this healthy young volunteer or a reconstruction artifact. This fact highlights the advantage of a complete process that starts from acquired image processing for 3D reconstruction of the vasculature of interest to reconstruction of virtual images. Then, by comparison between acquired and reconstructed images, one can assess the cause of any constriction along the vascular circuit and hence eliminate or confirm an artifact.

Newtonian assumption

A Newtonian model was adopted for blood rheology. In addition to the arguments developed in Section 3, blood flow simulations in cerebral aneurysms were not sensitive to blood rheology models [START_REF] Gambaruto | Sensitivity of haemodynamics in patient specific cerebral aneurysms to vascular geometry and blood rheology[END_REF]). Moreover, the blood viscosity model influences to a lesser extent blood flow computations than the modeling assumptions on the boundary conditions, especially inlet conditions, except when a precise computation of the wall shear stress is needed [START_REF] Chabannes | Hemodynamic simulations in the cerebral venous network: A study on the influence of different modeling assumptions[END_REF]).

Rigid walls

A relevant model describing the deformability of the cerebral venous wall remains a challenging issue for several reasons: (i) in vivo assessments of material properties are not perfectly handled; (ii) midsize vessels are not accessible to ultra-sound elastography for estimation of rheological parameters;

(iii) in vitro measurements in samples collected from dead individuals are accurately managed, but not performed in perfused living tissues, and (iv) the state of the art in fluid-structure interaction does not take into account the multilayered wall made up of composite materials. Therefore, even in recent 3D simulations of cerebrovenous flows such as those in jugular veins, walls are also assumed to be rigid [START_REF] Caiazzo | Computational haemodynamics in stenotic internal jugular veins[END_REF]).

Boundary conditions

This contribution is a first step towards multiscale models, in which 3D simulations are carried out in the region of interest and are coupled with reduced models of other compartments of the cardiovascular system, including the cardiac pump. To the best of our knowledge, only very few works are available in the case of the venous network (3D simulations of blood flow in the jugular veins coupled with a 1D model in major intracranial veins [START_REF] Caiazzo | Computational haemodynamics in stenotic internal jugular veins[END_REF]) and coupling of a 3D model with a simplified three-elements Windkessel model [START_REF] Chabannes | Hemodynamic simulations in the cerebral venous network: A study on the influence of different modeling assumptions[END_REF]), which represent in fact partial modellings. Moreover, the mechanical parameters are, most often, not obtained in the same individuals.

The extracranial compartment of the cerebral venous drainage is more deformable than the intracranial one [START_REF] Stoquart-Elsankari | A phasecontrast MRI study of physiologic cerebral venous flow[END_REF]). Cervical veins such as the upstream segment of the internal jugular vein that runs vertically down the neck side close to the internal carotid artery undergo significant deformations due to pressure pulses propagating through the carotid artery. These transient strains are here neglected because the present study focuses on the intracranial compartment, the virtual exit veins serving only for proper application of outlet boundary conditions. Artificial inflow and outflow boundaries are introduced by cutting the computational domain from the rest of the circulatory loop. The inflow velocity of a venous circuit has a quasi-steady profile, a sensitivity analysis of its magnitude being performed in [START_REF] Chabannes | Hemodynamic simulations in the cerebral venous network: A study on the influence of different modeling assumptions[END_REF]). The results showed that for cerebral veins, setting the inlet boundary condition has an important impact on the forces created by blood flow, thus highlighting the importance of deriving values for these conditions from clinically measured data at some probe stations in order to enhance the accuracy of the computed hemodynamical quantities of interest. However, technical limitations do not allow these measurements in small-bore vessels such as entry veins of the present model.

Simple coupling of the downstream vasculature of the neck and cranial part of the thorax using reduced models relies not only on electrical analogues (0D models), which yet should be used with proper estimation of personalised parameter values, but also adequate 1D model of deformable (both distensible and collapsible) veins to avoid reflection artifacts of pressure waves. Hence, a proper personalised vein law must be first determined. Furthermore, the unsteady venous outflow depends on the cyclic changes in internal and external pressures associated with the venous return in the right atrium of the beating heart under the influence of breathing. The body's position yields another constraint, as it affects the outlet drainage path and hence flow distribution among all possible exit vessels [START_REF] Gisolf | Human cerebral venous outflow pathway depends on posture and central venous pressure[END_REF]). Therefore, the present computational domain focuses on the intracranial compartment.

Computational domain

The computational domain must derive from a relevant geometry, in addition to an accurate topology. Slight errors in the geometric modeling can lead to erroneous blood flow distribution (Figure 4), the global computed pressure field being affected by a local geometry alteration that can result from anatomical influences or reconstruction artefacts. Comparison of original images to reconstructed ones will allow us to analyze the cause of any constriction along the vascular circuit and hence eliminate or confirm an artifact.

Concluding remarks

A methodology was presented for subject-specific blood flow simulations in the intracranial venous network to prove the feasibility of the pipeline from medical image acquisition to simulation. A new algorithm is currently under development and it will be used to carry out high-performance computing that allows the use of refined meshes. In addition, both anatomical and functional imaging (i.e., associated with MR velocimetry) is beginning to investigate 3D blood flow simulation in macrovascular networks.

One of the main open questions is linked to the validation of the results using phantoms or, better, in vivo measurements by functional imaging. However, measurement themselves have limitations due to space and time resolution of devices and vessel accessibility to measurements. The development of virtual angiography, which simulates the physical image acquisition sequences to create realistic synthetic data, provides a promising perspective that helps to interpret some features of the derived computational meshes [START_REF] Fortin | Extension of an MRI simulator software for phase contrast angiography experiments[END_REF][START_REF] Ford | Virtual angiography for visualization and validation of computational models of aneurysm hemodynamics[END_REF]; [START_REF] Jurczuk | A new approach in combined modeling of MRI and blood flow: A preliminary study[END_REF]).

Figure 1 .

 1 Figure 1. Cerebral venous network. 1: virtual exit veins; 2: superior sagittal sinus; 3: straight sinus; 4: confluence of sinuses; 5: lateral sinus (transverse portion); 6: lateral sinus (sigmoid portion); 7: vein of Galen; 8: internal cerebral vein; 9: basal vein; 10: superior cerebral veins; 11: superior anastomotic veins.

Figure 2 .

 2 Figure 2. Unsteady Navier-Stokes simulation (iteration 135 over 200, i.e. 0, 35 second in the second cardiac cycle): on the top N.a), the streamlines of the whole cerebral venous network, below on the left N.b), the velocity field on the confluence of sinuses and on the right N.c), the velocity field on the superior sagittal sinus (vectors represent the velocity direction and magnitude and colors represent the Y velocity component). Each Nth set a), b) and c) corresponds to the Nth mesh.

Figure 3 .

 3 Figure 3. Comparison between MRI measurements of flow rates (ml/s)(average and standard deviation) found in literature and our model in different locations of the intracranial venous network. SSS : Superior Sagittal Sinus Fig. 1 (2), St. S : Straight Sinus Fig. 1 (3), L. TS: Left Transverse Sinus Fig. 1 (5), R. TS: Right Transverse Sinus Fig. 1 (5)

Figure 4 .

 4 Figure 4. Side effects of geometric modeling inaccuracies on blood flow simulation relevance. A local vessel diameter reduction can lead to an abnormally elevated pressure field (normalized values maximum in red).

Table 1 .

 1 Venous geometry: approximate magnitude (mm) of the radius at some typical sites.

	Vessel name	Radius (in the literature)	Radius (MRA, present study)
	Right internal jugular vein	3.0 → 7.0	Ho et al. (2013)	2.59
		2.52 → 3.99 Müller, Toro (2014a)
	Left internal jugular vein	3.0 → 7.0	Ho et al. (2013)	2.34
		3.99 → 6.18 Müller, Toro (2014a)
	Superior sagittal sinus	2.3 → 3.0	Ho et al. (2013)	2.1 → 4.9
		2.0 → 3.19	Müller, Toro (2014a)
	Straight sinus	1.8	Ho et al. (2013)	3.3
		2.5	Müller, Toro (2014a)
	Right transverse lateral sinus 4 → 5	Schaller (2004)	2.9
		1.78 → 2.52 Müller, Toro (2014a)
	Left transverse lateral sinus	4 → 5	Schaller (2004)	4.9
		3.09 → 4.37 Müller, Toro (2014a)
	Right sigmoid sinus	2.52	Müller, Toro (2014a) 3.8
	Left sigmoid sinus	3.57 → 3.99 Müller, Toro (2014a) 5.2
	Vein of Galen	3.09 → 4.0	Müller, Toro (2014a) 2.06 → 2.50
	Internal cerebral vein	1.26	Müller, Toro (2014a) 2.07
	Basal vein	1.26	Müller, Toro (2014a) 2.03
	Superior cerebral veins	-		1.55 → 1.83
	Superior anastomotic veins	1.0	Müller, Toro (2014a) 1.07 → 1.86

Table 2 .

 2 However, patient-specific blood flow data are usually not collected from routine

	Sto 2 Re	and Sto = L * ω ν	are the Reynolds number, Strouhal
	number and Stokes number, respectively.

Data used in the present model to compute the flow-governing dimensionless parameters are given in

Table 2 .

 2 Values of the mechanical parameters in the model.

	Parameter	Value Unit	Source
	Mass density	ρ	1055	[ kg m 3 ]	Thiriet (2008)
	Dynamic viscosity	µ	3.5 •10 -3	[P a • s] Thiriet (2008)
	Kinematic viscosity	ν	(2.8-3.8) •10 -6	[ m 2 s ]	Thiriet (2008)
	Flow pulsation	ω	1	[ 1 s ]	Thiriet (2008)
	Cross-sectional velocity (clinical studies, jugular vein)	V *	(10-11) • 10 -2 (8.5-11.3) • 10 -2	[ m s ]	Ho et al. (2013) Ogoh et al. (2011)
			(30-50) • 10 -2		Schaller (2004)
	Cross-sectional velocity (clinical studies, superior sagittal sinus)	V *	15 • 10 -2	[ m s ]	Schaller (2004)
			(15.2 ± 3) • 10 -2	[ m s ]	Gideon et al. (1996)

Table 3 .

 3 Values of the dimensionless governing parameters. Hence, scarce literature data are used for the velocity magnitude in cerebral veins. Computed values for the Reynolds, Stokes, and Strouhal number are reported in Table3. The order of magnitude of the Reynolds number shows the importance of the convective forces, whereas the values of the Stokes and Strouhal numbers are important to assess the flow unsteadiness. The flow in the cerebral veins is thus governed by the Navier-Stokes equations for an incompressible viscous fluid, in a quasi-steady regime.

	Vessel name	Reynolds number Stokes number Strouhal number
	Internal jugular vein	75 → 1055	1.38 → 3.84	0.014 → 0.030
	Superior sagittal sinus 90 → 144	1.10 → 1.75	0.013 → 0.021
	clinical examinations.

Table 5 .

 5 Finally, pressure and velocity

	October

Table 5 .

 5 CPU Time for the different meshes

	Mesh 1 Mesh 2 Mesh 3 Mesh 4 Mesh 5
	CPU time (s) 3,447	4,853	8,386	13,747	43,122

www.vivabrain.fr

October 21, 2016 a journal

http://www.freefem.org

http://www.vtk.org

http://www.paraview.org


O Meste

M Meo

S Signore

M Rota

A new tool for the action potential repolarization dynamic analysis: Application to the discrimination of diabetic and control cells

In CinC2015 we have shown that the dynamics of the action potential (AP) repolarization could be tracked throughout the stimulation course. Despite some valuable outcomes, the populations of interest (control, Ctrl, and streptozotocin-induced, STZ, diabetic mice) could not be significantly distinguished in term of dynamics because of the global extracted feature. In this study, the computation of new features for each repolarization percentage allows an accurate and meaningful characterization of the two groups leading to a significant classification. APs in isolated left ventricular cardiomyocytes obtained from Ctrl and STZ mice were measured by patch-clamp. The progressive changes in AP repolarization for individual cells were tested on a set of 100 consecutive excitations at 2Hz pacing rate. The corresponding repolarizations are stacked in a matrix decomposed with a new approach. Observations are modeled as a sum of vectors multiplied by specific polynomial functions. This approach is similar to the Singular Value Decomposition (SVD), but the corresponding scalars are replaced by these functions. Model unknowns are estimated by using an alternated least square algorithm. Finally, the mean of the polynomial first derivative is computed for each repolarization percentage as a representative feature. A Wilcoxon signed rank test (p<0.05) has been applied on the features from the two groups. We can observe a significant difference in the late repolarization phase (70%-95% repolarization), with a singular behavior in correspondence with the AP profile shoulder onset (80%).

Introduction

The duration of the action potential (AP) in cardiomyocytes is an important variable controlling the electrical properties of the normal and pathological myocardium. Prolongation of the AP in myocytes may represent the ba-sis for the increased risk of arrhythmia with diabetes [START_REF] Tomaselli | Sudden cardiac death in heart failure: The role of abnormal repolarisation[END_REF], but the determinants of these abnormalities remain to be elucidated. We have shown in [START_REF] Meo | Reduction in Kv Current Enhances the Temporal Dispersion of the Action Potential in Diabetic Myocytes: Insights from a Novel Repolarization Algorithm[END_REF] that the diabetic condition is associated with alterations in the temporal dynamics of the AP profile in myocytes, a factor that may originate electrical instability. This assessment has been performed by using statistics computed over APs associated to isolated left ventricular myocytes obtained from control (Ctrl)and streptozotocin-induced (STZ) diabetic mice. In [START_REF] Meste | Diabetes affects the temporal dynamics of the repolarization properties of cardiomyocites[END_REF], these statistical properties (mean and variance) have been replaced by the analysis of the temporal dynamics of the repolarization phase. It has been shown that the two groups could be partly distinguished by introducing a new automatic analysis based on a transformation of each AP.

The objective of this paper is not restricted to only fully distinguish the two groups but also to produce additional predictive information on the ionic currents involved in the observed differences. For instance, if the groups are significantly different at specific percentages of repolarization, it implies that ionic currents active in these membrane potential ranges are affected by the hyperglycemic condition. To achieve this task, we introduced a new decomposition method, well adapted to the observation, based on the aforementioned AP transformation. It is shown that the proposed approach succeeded in distinguishing the electrical behavior of the two groups of cells, corroborating the impact of diabetes on repolarizing Kv currents.

Material

A group of female C57Bl/6 mice between 3 and 4.4 months of age was studied in accordance with the Guide for Care and Use of Laboratory Animals; experiments were approved by the local animal care committee (IACUC). Hyperglycemia was induced by injecting animals with streptozotocin (Sigma, concentration 5 mg/L) for approximately one week. Blood glucose levels were measured 6-10 days after the last injection of STZ. Mice with blood glucose concentration higher than 400 mg/dl were ascribed to the diabetic group. Furthermore, we examined only animals at less than 1 month after the onset of hyperglycemia, so as to minimize the interference of other factors (e.g., structural or mechanical changes) on myocardial electrical activity.

A set of 41 cells was isolated from left ventricle in 14 control mice, whereas the STZ group consisted of 76 cells obtained from 10 animals. Isolated cardiomyocytes were placed in a bath of Tyrode solution on the stage of an IX53 (Olympus) microscope for whole-cell patch-clamp measurements. APs were assessed in current-clamped myocytes through an Axoclamp 900A amplifier (Molecular Devices). Electrical signals were digitized using a 250 kHz 16-bit resolution A/D converter (Digidata 1550, Molecular Devices) and recorded through pCLAMP 9.0 10 software (Molecular Devices) with low-pass filtering at 2 kHz. Recording pipettes were fabricated using a horizontal (P-1000, Sutter Instrument) glass microelectrode puller. Current pulses 1.5 times threshold were applied for cell stimulation. A sequence of 100 APs has been triggered by the stimulation protocol at a 2 Hz pacing rate.

Method

The goal of the presented method is not to measure an absolute trend in the observation, for each beat and for each repolarization magnitude, but to determine how AP profile dynamics is globally affected throughout the stimulation sequence. Indeed, it can be seen as an extension of the singulat value decomposition (SVD) approach where an observation can be represented as the weighted sum of vectors, and each of its weights is not a single scalar value, but a function. For the SVD, the decomposition of each AP observation vector x i (i = 1, . . . , I) is defined by:

x i = a 1,i v 1 + a 2,i v 2 + a 1,i v 3 + . . . (1)
where the vectors v j are constrained to be orthogonal and normalized. This decomposition does not take into account any temporal structures hidden in vectors x i , meaning that if their components are shuffled the corresponding v j will be shuffled accordingly and producing the same a j,i .

The new model we propose consists in decomposing the observations such that :

x i = p 1,i • v 1 + p 2,i • v 2 + p 3,i • v 3 + . . . (2)
The symbol '•' stands for the Hadamard product corresponding to the element-wise product. The application of this decomposition to the AP matrix consists in considering a x i as the successive time values throughout the stimulations, for a given amplitude level. In this case we assume a temporal structure captured by both the p j,i and v j that will allow a lower complexity of the decomposition compared to (1). When using SVD, the variation information is often shared by different components, thus making its interpretation more difficult. In other words, it is expected that the new modeling would facilitate this task because of the lower number of summation terms. For our practical case, the model will be approximated by an order 1 decomposition:

x i = p i • v + e i (3)
with e i the approximation error. Assuming that the vector p i is fully parameterized by using a linear combination of basis vectors, that is:

p i (n) = K-1 k=0 b k (n)θ i,k (4)
with b k (n) a polynomial basis for instance, the minimization of

I i=1 e i 2
2 provides an estimate of the parameters of p i and v, with v 2 2 = 1. The latter constraint avoids indeterminacies of the scaling and imposes the uniqueness of the decomposition. This estimation is performed by using an alternated least squares approach with the two following solutions:

v = ( I i=1 diag(p i ) 2 ) -1 ( I i=1 (x i • p i )) (5)
and

p i = M v θi (6) θi = (M T v M v ) -1 M T v x i (7)
with M v = B • (v1I T ) and where the b k (n) are stored column-wise in B and 1I stands for the unit vector. Note that for the full extension in (2), by using similar approach to the SVD computation, all the unknowns can be estimated with the constraints that the summation terms are orthogonal (but not the v j ) and v j 2 2 = 1. Once the p i are estimated the mean p is computed as p = 1 I i p i and ṽ = p • v such that 1 I i x i = ṽ, neglecting the approximation error. This also allows to write x i = pi • ṽ with pi (n) = p i (n)/p i (n). This normalization step is crucial because the algorithm could converge toward -p i and -v instead of p i and v while the Hadamard product of the two vectors produces the same result. This is important in order to compare the derivatives of the pi (n) for different observation sets regardless this indeterminacy.

In our application each p i is modeled as an order 3 polynomial function with parameters to be estimated. The originality of the approach is to not assess the discriminative power of the ṽ, however distinct for each cell, but to focus on the property of the transformed polynoms pi . If for a given AP magnitude the polynom tends to decrease the trend of the time values throughout the time, then it should be observed through the mean of its derivative, regardless the profile of the corresponding ṽ. The main property of ṽ is to summarize all the common information shared by all the x i in a different manner than with the SVD. For a given cell and for each i value, according to each AP magnitude, the mean of the derivative of pi is computed and stored. The discrimination of the diabetic and the control cells will be based on these last values.

It has been shown in [START_REF] Meo | Reduction in Kv Current Enhances the Temporal Dispersion of the Action Potential in Diabetic Myocytes: Insights from a Novel Repolarization Algorithm[END_REF] that, according to the mean AP for the two groups, it is more relevant to normalize the observations to avoid bias for the groups discrimination. To achieve this task, the previously computed parameters are normalized accordingly.

Results

The idea is to distinguish the control (41 records) from the diabetic (76 records) myocells specifically on the repolarization phase. First, for each AP a transformation is applied in order to get a strictly decreasing curve, which in turn is also transformed by computing its inverse [START_REF] Meste | Diabetes affects the temporal dynamics of the repolarization properties of cardiomyocites[END_REF]. Then for each stimulation sequence (the 30 first have been discarded) we get a matrix of 70 inverse repolarizations (see a sequence of stimulated control cell in Figure 1 and the corresponding 3D visualization of the corresponding matrix in Figure 2). Note that the horizontal axis stands for the normalized amplitude and the vertical one for the time. The application of the proposed decomposition to the APs consists in considering a x i as a column of the matrix.

The derivatives, computed over the stimulation index, of the corresponding pi , are given in Figure 3 and will be used to characterize the trend as previously explained. A ranksum test has been computed over the mean of the derivatives in order to discriminate the two groups and is given in Figure 4. In order to interpret the intervals of significance of this test, the mean AP and its corresponding derivative is plotted for the two groups in Figure 5. Finally, the median value of the discriminative parameter is computed for the two groups in Figure 6.

Discussion and conclusion

In previous works [START_REF] Meo | Reduction in Kv Current Enhances the Temporal Dispersion of the Action Potential in Diabetic Myocytes: Insights from a Novel Repolarization Algorithm[END_REF], [START_REF] Meste | Diabetes affects the temporal dynamics of the repolarization properties of cardiomyocites[END_REF] applied to a subset of the records it has been shown that automatic AP analysis can be performed accurately. It has not been demonstrated that the dynamics of the repolarization phase throughout the stimulation allows to distinguish the two populations while some specificity has been found. This could be due to the inherent random variability of the repolarization regardless the groups. While the classification of the two groups is not really relevant in this context, the related tools could reveal new insights of the diabetes effects over the repolarization. From the example given in Figure 3, it is clear that the multiplicative polynoms p i convey complex infor- It is noticeable that the two groups are significantly different only around the shoulder of the mean AP (see Figure 5 for localization) and mostly in the late repolarization phase. Referring to [START_REF] Nerbonne | Molecular Physiology of cardiac repolarization[END_REF], this localization could be due to changes of the Kv currents profile. In addition to the work in [START_REF] Meo | Reduction in Kv Current Enhances the Temporal Dispersion of the Action Potential in Diabetic Myocytes: Insights from a Novel Repolarization Algorithm[END_REF], where it is shown that these changes are due to the effect of diabetes, this new findings could help to identify which Kv currents are the most affected. Furthermore, the median value shown in Figure 6, brings additional information on the dynamics. The diabetic group exhibits a shortening trend (negative values) mostly for the late repolarization phase and conversely for the early phase. For the control group a shortening is observed only in the interval [80-65]% corresponding to the location of the mentioned shoulder. This result extend the observations given in [START_REF] Nobe | Chronic diabetes mellitus prolongs action potential duration of rat ventricular muscles: circumstantial evidence for impaired Ca2+ channel[END_REF] where similar behaviors are shown to be frequency dependent.

In order to investigate the role of specific Kv currents involved in the observed electrical alterations, new targeted experiments are needed. Alternatively, in silico models

Figure 1 .Figure 2 .

 12 Figure 1. 70 superimposed repolarizations of a control cell

Figure 3 .

 3 Figure 3. The derivative with respect to stimulation index of the multiplicative estimated polynomials.

Figure 4 .

 4 Figure 4. Comparison of the two populations (Ranksum). Significant when ≤0.05.

Figure 5 .Figure 6 .

 56 Figure 5. Mean and derivative (magnified by 100) of the control and diabetic populations. Note the shoulder shape in the interval [80-70].

Address for correspondence:

Pr. Olivier MESTE Laboratoire I3S -CNRS -UNS 2000 route des Lucioles 06903 Sophia Antipolis cedex, FRANCE E-mail address: olivier.meste@unice.fr


Olivier Meste

Stef Zeemering

Joël Karel

Theo Lankveld

Ulrich Schotten

Harry Crijns

Ralf Peeters

Pietro Bonizzi

Noninvasive Recurrence Quantification Analysis Predicts Atrial Fibrillation Recurrence in Persistent Patients Undergoing Electrical Cardioversion

This study proposes a novel noninvasive ECG-based approach to accurately quantify differences in the atrial activity (AA) signals of atrial fibrillation (AF) patients characterized by very similar global AF organization.

High-density body surface potential maps were recorded in 63 patients in persistent AF prior to electrical cardioversion (32 recurrences occurred after 4-6 weeks). AF substrate complexity was quantified as the degree of overall temporal organization of the AA pattern, determined from multi-variable unthresholded recurrence plots (RP).

Results showed that non-recurrent patients were characterized by temporally more stable propagation patterns than recurrent patients, suggesting a lower complexity of the AF substrate in the former. This was used to discriminate between recurrent and non-recurrent AF patients following electrical cardioversion with an accuracy of 72%.

Noninvasive assessment of AF substrate complexity in patients characterized by very similar levels of AF progression was improved by applying RP-based analysis on extracted AA signals. The proposed approach may help improving AF treatment selection and success.

Introduction

ECG-based noninvasive quantification of atrial fibrillation (AF) substrate complexity is increasingly being used to discriminate among AF types and to predict therapy outcome at any stage of the disease [START_REF] Lankveld | The ECG as a tool to determine atrial fibrillation complexity[END_REF]. In particular, predicting acute and long-term success of AF treatments may improve treatment selection and success, and prevent future cardiovascular complications. Many ECG-based AF complexity parameters have been suggested over the last decade and some first steps towards a standardized approach have been taken [START_REF] Bonizzi | Systematic comparison of non-invasive measures for the assessment of atrial fibrillation complexity: a step forward towards standardization of atrial fibrillation electrogram analysis[END_REF]. A limitation of those parame-ters is that they generally quantify stationary AF substrate properties and appear to be suboptimal in describing subtle differences in patients characterized by very similar AF substrate complexity. However, patients showing similar AF substrate complexity still may respond differently (on a short and long term) to the same treatment. This asks for more accurate ECG-based noninvasive methods able to unveil differences in the global AF organization of those patients, which are more subtle and more difficult to quantify [START_REF] Uldry | Measures of spatiotemporal organization differentiate persistent from long-standing atrial fibrillation[END_REF], and likely due to structural strengthening of already existing remodelling processes.

In this study, we introduce a novel approach based on unthresholded recurrence plots (RP) to describe the oscillatory patterns of atrial activity (AA) signals during AF and well represent the continuous spectrum of AF progression. The proposed approach is also applied to noninvasively predict AF recurrence in persistent AF patients undergoing electrical cardioversion.

Methods

BSPM data and pre-processing

Body surface potential maps (BSPMs) were recorded in 63 patients in persistent AF using a custom configuration of 184 leads with 120 anterior and 64 posterior leads (Ac-tiveTwo BSM Panels Carbon Electrodes, Biosemi B.V., The Netherlands; Figure 1). All patients underwent electrical cardioversion, and 32 patients showed AF recurrence 4 to 6 weeks after procedure. ECGs were sampled at 2048Hz, and downsampled at 256Hz. A one-minute segment was selected for each subject, low-quality leads were excluded (low signal-to-noise ratio, poor electrode contact, motion artefacts), and Wilson's Central Terminal was subtracted in line with conventional ECG analysis. After band-pass filtering the signals between 1 and 100Hz (3rd Figure 1. BSPM electrode configuration, comprising 120 anterior leads and 64 posterior leads. Thick circles: default positions of the precordial leads V 1 , . . . , V 6 . order Chebyshev), QRST cancellation was performed using an adaptive singular value decomposition method, inspired by the approach in [START_REF] Alcaraz | Adaptive singular value cancellation of ventricular activity in single-lead atrial fibrillation electrocardiograms[END_REF], with multiple QRST window templates defined using hierarchical clustering. The extracted atrial signals were post-filtered with a zerophase notch filter at 50Hz to suppress power line noise, and with a 3Hz zero-phase highpass filter (3rd order Chebyshev) to remove low-frequency residuals not related to (persistent) AF. Figure 3

Multi-variable recurrence plot construction

RPs provide a way to visualize and quantify recurrent behavior of the phase space trajectory of a dynamical system [START_REF] Eckmann | Recurrence Plots of Dynamical Systems[END_REF]. RPs can be generated from trajectory data. In this study, given the matrix X of size `⇥ N collecting all extracted AA signals from a patient (`= 184 leads, and N is the number of samples), each column of X was assumed to provide a vector x(i), i = 1, . . . , N, which represents a point in the phase space of the system. This allows to process all AA signals extracted from a patient simultaneously, by means of a multi-variable RP.

Traditionally, the Euclidean distance is used to compute the proximity of two points in the phase space. However, we replaced this distance with a similarity measure, namely the cosine of the angle between the vectors associated with two points in the phase space:

d(x(i), x(j)) = x(i) T x(j) ||x(i)|| 2 ||x(j)|| 2 (1)
This provides a sort of sample correlation which is not affected by the magnitude of the vectors, and allows to investigate the temporal recurrence of global AF spatial patterns over the all body surface covered by the BSPM electrodes.

Recurrence is then expected to happen not when the dis- tance is zero, but when similarity is high. Finally, by computing the similarity (1) between each pair x(i), x(j), an unthresholded RP matrix of size N ⇥ N was obtained for each patient, which was used for further analyses.

Recurrence Quantification Analysis of multi-variable AA RPs

RQA uses RPs to investigate and quantify the recurrent behavior in the phase space trajectory of a dynamical system. Under the hypothesis that an AA signal during AF is characterized by a certain stationarity, the corresponding AA cardiac dipole can be assumed to pass through the same points (or close to) after a period T . At such moments in time the distance d(x(i), x(i + T )) (with T expressed in samples) defined in ( 1) is expected to be 1 (or close to 1). With this assumption of stationarity fulfilled, averaging over a diagonal of the RP would provide a normalized autocorrelation value at a lag given by the index of the diagonal (for instance, averaging over the main diagonal would provide a normalized autocorrelation at lag zero).

However, it is also sensible to assume that this stationarity is likely not to last over a long time interval, and to expect breaks in the phase of the AF patterns to occur (piecewise stationarity). The number of those breaks is expected to depend on how much progressed AF is, and on the overall status of a patient. The more progressed AF becomes the larger the number of breaks in the phase, and the more nonstationary the AF patterns are expected to be. There is no analytical way to compute the average length of a piecewise stationary window in a patients (the time interval between two break points in the phase of the AF patterns). The length of this window then becomes a parameter of the method. There is a tradeoff on the window length: the shorter the window is the better nonstationarities can be handled, but the larger the variance of the estimated normalized autocorrelations over different windows. On the contrary, the larger the window is the more accurate the estimates of the normalized autocorrelations are, but the less it can adapt to nonstationary behaviors. In this study, we empirically determined that a window length of 500 samples worked well in terms of discriminating between the recurrent and non-recurrent AF patients.

After the window length was chosen, each RP was split into a sequence of non-overlapping diagonal blocks, all of the same size, each block containing 500 diagonals of equal length (see an example in Figure 2(a) for three sequential blocks). That means that by averaging over all diagonals of one block, all lags between 0 and 499 samples were assessed. In this way, we could determine a new autocorrelation function for every block, and adapt to nonstationarities. A normalized autocorrelation function was then obtained for each block (see Figure 2(b)). For a window length of 500 samples, twenty-four blocks were obtained from the RP of each patient, thus providing twentyfour normalized autocorrelation functions per patient.

The envelope of each autocorrelation function was then computed as the amplitude of its Hilbert transform (thick solid lines in Figure 2(b)). This allows to work with smoother signals, not affected by discontinuities in the phase. Finally, for a given patient the median of all the twenty-four envelopes at each time lag was computed, providing a 500 sample long signal, r(p), with lags in the range p = 0, . . . , 499 samples (Figure 2(c), thick line overlapped to the three envelopes). This signal includes information about the recurrent behavior present in the original data, at the same time accounting for nonstationary properties of this behavior. Two analyses were then carried out based on the signal r(p).

3.

Analyses and Results

Qualitative analysis of AF nonstationary behavior

Figure 3(a) shows the r(p) signals for all 63 patients. A visual inspection of those curves shows that they are characterized in general by two distinct behaviors: an early phase characterized by a decreasing autocorrelation value (for low values of the time lag p), and a later phase characterized by an approximately constant behavior, with very low autocorrelation values (for larger values of the lag p). To highlight this common behavior and partly remove intra-individual differences, each r(p) signal was normalized by the corresponding median value of the constant phase (calculated as the median of r(p) in the interval 150  p  450; Figure 3(b)). Finally, given the two groups composed by recurrent and non-recurrent AF patients, respectively, the mean of all r(p) signals was computed for each group (over each time lag; Figure 3(c)). A visual inspection of the two resulting curves shows that both groups are characterized by the two distinct behaviors highlighted above. Moreover, it can be noticed that the mean curve of the non-recurrent group shows values that are larger than those of the recurrent curve at the early phase (10  p  60). That means that in this early phase, the non-recurrent patients show in average a larger autocorrelation, which can be interpreted as a higher regularity of the corresponding global AF patterns, which tend to pass for longer from the same locations. This is in line with the expectation that electrical cardioversion is more likely to be successful in those patients with less electrostructural remodeling of the myocardium, reflected in a lower complexity of the AF substrate and more organized AA patterns.

3.2.

RP-based classification of recurrent and non-recurrent AF patients

Based on the previous qualitative analysis, an open question is: what is the value of the time lag p which allows to separate the recurrent and non-recurrent patients best? To answer this question we first computed the cumulative function of each normalized signal r(p) (Figure 4). This allows, given a specific time lag, to take into account the contribution of the previous lags as well. We then performed a leave-one-out cross validation on each time lag to determine the most discriminant lag. We found that 200 samples was the value of the time lag which gave the highest accuracy in terms of classification of recurrent and non-recurrent AF patients (72%; AUC 0.58). When looking back at Figure 3(c), it can be noticed that p = 200 is indeed at the end of the early decreasing phase, before the curves become constant, underlining that all discriminant power is exploited by this value.

AF complexity assessment

For completeness, AF complexity was also assessed on the extracted AA signals by means of fourteen state of the art noninvasive time and frequency domain parameters that appeared frequently in the last decade of non-invasive AF complexity literature [START_REF] Bonizzi | Systematic comparison of non-invasive measures for the assessment of atrial fibrillation complexity: a step forward towards standardization of atrial fibrillation electrogram analysis[END_REF]. A generalized linear model was individually applied to each complexity parameter. None of those parameters showed to be a significant discriminant of AF recurrence.

Conclusions

This study proposed a novel RP-based noninvasive quantification of AF substrate complexity. Novelty resides in the use of a normalized autocorrelation instead of the Euclidian distance to build RPs and in extracting from the RP matrix a signal that contains information about recurrences in the original data, accounting for possible nonstationaries of those recurrences. In particular, the proposed approach showed to be able to unveil more subtle information about AA patterns in persistent AF patients, which could be used to discriminate between recurrent and non-recurrent AF patients following electrical cardioversion with an accuracy of 72%. Non-recurrent patients were characterized by temporally more stable propagation patterns than recurrent patients, suggesting more organized AA patterns and a lower complexity of the AF substrate in the former, in line with what could be expected from a pathophysiogical view-point. On the contrary, state of the art noninvasive parameters for the assessment of the complexity of the AF substrate did not show to be able to discriminate between recurrent and non-recurrent AF patients, suggesting that more sophisticated tools are needed to describe patients characterized by very similar levels of AF organization. Main findings of this study include that: persistent AF is a complex nonstationary process, characterized by different nonstationary behaviors on the short and long term; following electrical cardioversion of persistent AF, non-recurrent and recurrent patients show different nonstationary behaviors, and non-recurrent patients are characterized by more organized AA propagation patterns; differences in the nonstationary behaviors of those two groups can be described by the proposed approach and used to obtain a good classification. Finally, the method proposed in this study provides a more accurate description of the AA oscillatory patterns hidden in AF ECG recordings and helps predicting AF recurrence in persistent patients undergoing electrical cardioversion.

  (a) shows the r(p) signals for all 63 patients.

Figure 2 .

 2 Figure 2. The main phases of the proposed approach. Read Section 2.3 for details.

Figure 3 .

 3 Figure 3. The three steps of the analysis of AF nonstationary behavior. 3(a), raw median r(p) signals from all 63 patients. 3(b), the normalized r(p) signals. 3(c), the mean signals for the non-recurrent and recurrent AF groups.

Figure 4 .

 4 Figure 4. Cumulative functions of the normalized r(p) signals in Figure 3(b).


Samuel Huron
email: samuel.huron@cybunk.com

Écom

Sheelagh Carpendale
email: sheelagh@ucalgary.ca

Jeremy Boy

Jean Daniel Fekete
email: jean-daniel.fekete@inria.fr

Using VisKit: A Manual for Running a Constructive Visualization Workshop

Figure 1: The procedure of VisKit: 1/ Participants are instructed by the workshop crew; 2/ Groups of participants encode data with tangible tiles; 3/ Participants fill a questionnaire to assess their visualization; 4/ Participants change their visualizations for communication purposes; 5/ Participants read, critique, and assess each other's visualizations; and 6/ The workshop crew conclude the workshop with discussion.

INTRODUCTION

As one response to the continuing spread of infovis online [START_REF] Viegas | Manyeyes: a site for visualization at internet scale[END_REF], in newspapers [START_REF] Segel | Narrative visualization: Telling stories with data[END_REF], and on TV [START_REF] Huron | Bubble-TV: Live Visual Feedback for Social TV Broadcast[END_REF], we have been running workshops to teach infovis to the general public. An outcome of these workshops is VizKit, a toolkit designed for teaching infovis to a board spectrum of people in a two hour workshop. In this paper we describe VizKit and the procedure we use during these public workshops. While running these workshops only reaches a few people each time, the people who attend have been enthusiastic and have, for instance, chosen to stay for the full length of the workshop rather than take part in a public fair. While it is well recognized that teaching infovis in a classroom setting is challenging [START_REF] Hearst | Vis, the next generation: Teaching across the researcher-practitioner gap[END_REF][START_REF] Kerren | Teaching information visualization[END_REF], if we want to democratize visualization to lay audiences it is also important to think about teaching outside classrooms [START_REF] Bhargava | Designing Tools and Activities for Data Literacy Learners[END_REF][START_REF] Duarte | Handmade visual[END_REF].

Both the toolkit, VizKit, and the procedure are based on constructive visualization principles [START_REF] Huron | Constructive visualization[END_REF]. We used these principles be-cause in theory they make learning more accessible [START_REF] Huron | Constructive visualization[END_REF], and when studied [START_REF] Huron | Constructing visual representations: Investigating the use of tangible tokens[END_REF] it appears that this might be the case. Since our intention was to provide an inclusive, accessible, and reflexive experience for lay audiences, when designing our authoring tools, we focused on supporting dynamic and collaborative visualization authoring. We worked towards minimizing prerequisite skills; avoiding technological restrictions; bringing the cognitive processes to the fore; and leveraging peer learning -an approach in which people can learn from each other.

In this paper we describe the VizKit components and outline how we have been using VisKit to run two hour constructive visualization workshops. Our procedure was specially designed to provide a first experience of data mapping, including both authoring and interpreting phases (See Fig. 1).

THE VIZKIT COMPONENTS

We designed this workshop for 15 people. We find that it works best if we can have a 2 person crew. While a projector and laptop are only used for a short period of time, it is definitely useful to have them. Then we have some prepared printed material and the physical authoring toolkits. Fig. 2 shows VizKit's components.

Printed Material

Instructions. There are three types of instruction sheets: 1) five datasets to visualize; 2) five personas (fictional characters) that describe the owners of each dataset and a set of questions about the features expressed in each dataset; and 3) a peer assessment grid to help people when they are reading the visualization. Datasets. We selected a type data that most adults can understandcategorized bank expenses. The datasets present 5 categories over 4 months. The categories are: Entertainment, Bar and Restaurants, Groceries, Transportation, Travel, and Rent. Each value is rounded to a multiple of $25 to simplify calculations. The different trends and features in the data such as increases, decreases, outliers, voids, and/or stability have been carefully crafted to correspond to the expenses of individual personas. We have five unique datasets, where each dataset has its own recognizable data features. Persona Sheets and Self-Evaluation Questions. These sheets provide a short story about each persona. The story explains the trends and/or features in the data and provides some cultural context, making it more likely that participants will "relate" to the data. Following the story is a set of questions (5 or 6 on average, each having two parts) about the different trends and/or features in the data. These questions ask participants whether they have identified a given feature, and whether they have made this feature visible in their visualization (dichotomous answers, i. e., yes/no). A Peer Assessment Grid. This is a table with questions as rows and participants as columns. The questions are the same as those in the self-evaluation sheet mentioned above. They relate to the same features for each dataset. The rationale behind the design of this sheet is to ease into a peer assessment activity in which participants try to read the visual representations created of other participants. Thus the participants engage in try to identify data features in visualizations authored by other the participants.

Authoring Toolkit

We chose to use tangible tiles similar to those used by Huron et al. [START_REF] Huron | Constructing visual representations: Investigating the use of tangible tokens[END_REF] since it was noted with these tiles lay audiences readily authored visual representations, dynamically and expressively. Therefore the workshop author toolkit is based on tokens. It includes two boxes, A and B, that each contain four cases that hold different colored tiles (Fig. 2). Each case contains thirty-six plastic tokens-25 * 25 * 5 millimeters in volume. Box A contains orange, purple, yellow, and black tokens; and box B contains red, green, blue, and white tokens (see Fig. 2). For instructors who wish to create their own toolkit, the 10 sets of boxes and tokens took three crew members and about two days of work using a laser cutter, glue, and hammer. We released the blueprint and documentation for the kit as open hardware 1 . Ideally you will have a space large enough to hold 6 large tables at which groups of 3 people can be comfortably seated. Before starting the workshop, setup your 5 small group workspaces with a table and 3 chairs. Also set up the crews workspace on a separate table. We also recommend setting up one table per group. You can distribute one toolset by table, and one dataset by table. Keep the other instruction sheets to distribute at the right moments. Labelling each table with a group number is also useful. If possible use a circular layout for the tables, as this will simplify the circulation of people during the critiquing phase of the workshop (see Fig. 3).

RUNNING THE WORKSHOP

The workshop consists of six phases: 1) introduction, 2) authoring, 3) reading and in group peer-evaluation, 4) editing for communication, 5) critiquing and cross group peer-evaluation, and 6) conclusion and wrap-up (see Fig. 1). This phase consists of a presentation preamble that introduces the components and procedure of the workshop (see Fig. 4). The participants are asked to form groups, taking seats around tables on which a dataset and the authoring kits have been placed. The datasets are described as bank statements belonging to an imaginary friend. We introduce the task to the participants by asking them to help their imaginary friend, who is having trouble with expenses. Note that the actual personas -that is the details of what the troubles areare not introduced in this phase. Each of groups is given a dataset, Fig. 6, and invited to use the tiles to create a visualization to help the fictional friend to manage their expenses. To simplify the cognitive processing, we suggest that one tile correspond to a base unit of 25. The intent of this authoring phase is to allow people to actively think about how to transform data into a visual form through practice (see Fig. 5). The dataset, the authoring tool, and the task are kept simple to encourage the group to create a visual representation. This process often triggers questions in the group on what and how to represent the data. In Phase 3, the persona sheets containing questions are handed out. Participants are asked to read it and reply to the questions. Participants are asked to relate the data to the fictional character they are supposed to be helping out (Fig. 7). The self-evaluation questions are designed to help them understand the trends and/or features in the data, and allow them to establish which of these they have already identified and represented in their visualization(s), and those they have not. The design rationale of this phase is to introduce a first external constraint: Does the visual mapping represent the main data features? The goal is to prompt people to take a reflexive moment about their authored visualization. To first re-read the visual representation and the dataset at the light of the persona. Then to identify main data features in the dataset. Finally they discuss collaboratively whether or not they have represented these features in their visualization. Here is a sample persona.

Phase 1: Introduction

Phase 2: Authoring

Phase 3: Reading and in Group Evaluation

The student is not originally from this city, but he has come here to study. Every year, he moves out of his apartment for the summer break, and goes back home to where his parents live. There, he spends very little money on groceries, and a lot on recreation. At the beginning of the Fall term, he moves back to school, finds a new apartment, and pays a new deposit. He also needs to pay for his transport card, which will be sufficient for all commuting during the year. He starts buying groceries again, and spends less on recreation. Finally, he spends quite a lot of time going out, and spends a relatively large amount of his money there.

Here are the associated sample questions for this persona: In Phase 4, participants are asked to reconsider their visual representation(s) for purpose of communication, i. e., as a "stand-alone" piece that the persona should be able to understand without their assistance. They are also invited to edit their visual representation(s) in response to the information provided on the persona sheet and self-evaluation questions (Fig. 8). Note that for this phase, participants are allowed to use other stationary equipment like pens, paper, and post-its. This phase introduces two other challenges: How to transform the visual representation, to make it readable for someone else? and How to correct the visual representation according to the persona sheet? The rationale behind this phase is to encourage participants to correct the visual representation, to reflect on the data features, and also to think about how to adjust their representation so that it will be more communicative.

Phase 5: Critiquing: Reading Each Other's Visualizations

In Phase 5, the peer-evaluation sheet is handed out, and participants are asked to nominate an evaluator for their group, i. e., a groupmember who stays seated next to the group's visualization; and one (or more) visitor(s). The visitors are instructed to go from one table to the next to try to interpret the different visualizations that have been created (see Fig. 9). They are asked to orally describe the different trends and/or features they can identify to the evaluator. The evaluators are instructed to keep track of the different trends and/or features listed in the peer-evaluation questionnaire that each visitor comes to identify, by checking the appropriate cell in the table (see Figure 10). Note that a evaluator is not allowed to discuss the visual encodings with a visitor from another group, until it is considered that the latter has identified all the insights shown by the visualization. Then they can spend some time explaining and critiquing the design and encoding choices that were made. In addition, participants are made aware of the dataset duplicates, i. e., that pairs of groups had actually been working with the same datasets. These duplicates are used to create a mini-challenge, in which participants are asked to identify which group worked on the same dataset they did. This helped participants to realize that there are a variety of encoding possibilities, and therefore that there is space for creativity offered by visualization design, even when a base unit is fixed.

Phase 6: Conclusion

In Phase 6, all personas and datasets are presented to the different groups, the different trends and/or features in each dataset are described using previously designed visualizations, and dataset duplicates are revealed (Fig. 11). To finish, some common visualization types are presented (e. g., scatterplots, node-link diagrams, etc.,) to emphasize the array of possibilities.

CONCLUSION

We have presented VisKit complete with detailed instructions on how to use VisKit to run a constructive visualization workshop. These workshops were designed to engage non-experts in a first experience of mapping data to a visual representation in less than 2 hours. VisKit lets people experiment with authoring, reading, editing, and critiquing a visualization. We have used VisKit in various environments and have received positive feedback. By presenting VizKit together with these instructions, which can be used as an instruction manual, we hope to engage more people in both exploring spreading visual literacy among the general public and the potential of using constructive visualization for learning purposes.

Figure 2 :

 2 Figure 2: The components needed for running the workshop.

1Figure 3 :

 3 Figure 3: Setting up the workshop.

Figure 4 :

 4 Figure 4: Presenting: Explaining how the workshop will run.

Figure 5 :

 5 Figure 5: The group author a visualization with the authoring kit .

Figure 6 :

 6 Figure 6: Example of a given dataset.

Figure 7 :

 7 Figure 7: Reading the visualizations the group has created.

1 . 4 . 4 4 :Figure 8 :

 14448 Figure 8: Editing: adjusting the visualization so that data features can be read more easily.

Figure 9 :

 9 Figure9: Critiquing: the person with the visitor role is trying to read the visualization of another group, while the person who is the evaluator is trying to detect if the visitor is discovering the encoded features.

Figure 10 :

 10 Figure 10: Example of a Peer assessement questionnaire for group1.

Figure 11 :

 11 Figure 11: Conclusion: wrapping up the workshop with group discussion

ACKNOWLEDGEMENTS

We want to thank all the people who helped with the work of creating VizKit, and designing and running workshops: Mathieu Legoc, Romain Di Vozzo, Sylvia Frederikson. We also thank the people who invited us to different events: Charles Nepote, Sylvia Frederikson.