Dominique Chapelle

Luca Formaggia

Philippe Helluy

Yvon Maday

Kent-Andre Mardal

Sébastien Martin

Christophe Prud'homme

Stéphanie Salmon

Christophe Prud'

Mathematical modeling, analysis and simulations for fluid mechanics and their relevance to in silico medicine

Keywords:

Mathematical modeling, analysis and simulations for fluid mechanics and their relevance to in silico medicine

Applications to the study of the coupled eye-cerebral system 2.1 Design and development of a lumped model for the coupled eye-cerebral system . . .

Préambule

Ce manuscrit rassemble mes contributions portant sur le développement de nouvelles méthodes mathématiques et numériques pour l'analyse des écoulements biologiques, considérés comme des phénomènes complexes multi-physiques et multi-échelles. La description des mécanismes sousjacents découle des principes de base de la mécanique des fluides et se traduit dans des systèmes d'équations aux dérivées partielles ou d'équations différentielles. L'objectif global de ce travail est l'étude de ces équations aux niveaux continu et discret et de leur couplage, ainsi que la mise en place d'un environnement de calcul permettant l'implémentation fiable et efficace des différentes méthodes numériques proposées pour approcher les solutions de ces problèmes. Les simulations numériques ainsi développées prennent en compte des géométries réalistes et une validation détaillée est proposée à travers des comparaisons avec des données expérimentales, en vue de leur application à des questions biomédicales spécifiques. L'approche in silico, consistant à combiner la modélisation basée sur des données et la modélisation fondée sur les lois de la physique, dans le même environnement de calcul, est devenue de plus en plus populaire dans le domaine des sciences de la vie ces dernières années et a motivé le développement de nouveaux concepts mathématiques et numériques. Malgré des avancées significatives dans la modélisation in silico de la physiologie humaine, l'étude de la dynamique complexe qui régit l'interaction entre différents fluides dans le corps humain suscite encore des questions extrêmement difficiles. La description des fluides biologiques fait intervenir une large gamme d'échelles spatiales et temporelles, entre le niveau moléculaire et celui des réseaux de quelques mètres, avec une durée qui peut aller d'une seconde pour le cycle cardiaque à plusieurs dizaines d'années pour une vie. La dynamique de ces fluides est influencée par l'interaction avec les tissus environnants, d'où la nécessité d'une approche multi-physique dans la modélisation. La représentation géométrique est très complexe et peu de données réalistes sont disponibles. En outre, dans les études expérimentales, de nombreux facteurs entrent en jeu et il est difficile d'isoler l'effet de chaque contribution sur les caractéristiques de l'écoulement. Dans cette direction, l'analyse statistique offre de nombreuses techniques pour aider à identifier des corrélations entre les données, mais la compréhension des relations de cause à effet à travers ces méthodes reste limitée. Dans cette perspective, la modélisation mathématique et numérique, basée sur les principes fondamentaux de la physique et adaptée à un cadre biologique, peut aider à identifier des relations de cause à effet dans l'action conjuguée de multiples facteurs et ainsi améliorer notre compréhension de ces phénomènes complexes ; bien évidemment, cette analyse est amenée à se développer en synergie avec d'autres approches.

Afin de mieux comprendre certaines difficultés rencontrées dans cette étude, il nous paraît important d'aborder à ce stade la question suivante : quelles sont les différences entre la modélisation des fluides biologiques et celle des fluides "classiques" ? Comparés à des fluides "classiques", ces fluides remplissent des fonctions biologiques essentielles (par exemple le transport de nutriments et 1 2 d'oxygène, l'interaction avec les cellules environnantes), qui doivent être prises en compte dans un modèle mathématique. Contrairement aux fluides "inertes", leurs propriétés changent au cours du temps, comme conséquence de différents processus biologiques liés à l'âge, à la maladie, etc. De nombreux défis sont engendrés par leur étude in vivo, car la plupart du temps, les fluides biologiques sont difficilement accessibles dans leur milieu naturel ; en outre, ils montrent souvent des propriétés et des comportements différents lorsqu'ils sont examinés in vitro ou in vivo. Par conséquent, l'extension des théories classiques de la mécanique des fluides et de l'approche mathématique correspondante à l'étude des systèmes biologiques soulève de nombreuses questions et représente actuellement un domaine de recherche très actif.

Tous ces aspects conduisent à : (i) de nombreuses difficultés liées à la dérivation d'un modèle mathématique pertinent et bien posé à partir de la littérature biologique et clinique ; (ii) un besoin important de développer de nouvelles approches théoriques et méthodes numériques ; (iii) une nécessité croissante de concevoir et de mettre en oeuvre un environnement de calcul flexible et efficace ; (iv) des efforts de recherche importants dédiés à la vérification, la validation et la prise en compte des incertitudes. Ces aspects sont cruciaux pour garantir que le modèle et ses solutions numériques sont significatifs du point de vue biomédical, permettant ainsi l'utilisation du cadre résultant comme un laboratoire virtuel où des idées peuvent être testées et de nouvelles hypothèses peuvent être formulées.

Dans ce contexte, les contributions décrites dans ce manuscrit se retrouvent à la croisée des mathématiques, de la modélisation et de la médecine. Comme résumé ci-dessous et détaillé dans les chapitres suivants, nous présenterons le développement de nouveaux modèles mathématiques et méthodes numériques, illustrés par des applications biomédicales. À travers les publications sélectionnées, chaque chapitre décrit les difficultés rencontrées et les résultats obtenus en réponse à différentes questions théoriques et pratiques et se termine par une discussion générale et des perspectives. Tous ces développements n'auraient pas été possibles sans une étroite collaboration avec des collègues mathématiciens, physiciens, informaticiens, biologistes et médecins. Leurs noms apparaissent dans la liste des références. Je mentionne en particulier dans cette introduction les contributions de plusieurs doctorants que j'ai pu encadrer, de manière formelle ou informelle et qui ont donné lieu à des publications communes.

Les chapitres sont structurés comme suit.

Dans le chapitre 1, nous nous intéressons aux modèles tridimensionnels, dans lesquels le mouvement d'un biofluide dans une géométrie complexe et réaliste est régi par les équations de Navier-Stokes sur un domaine faisant intervenir des conditions aux bords d'entrée et de sortie. Différents aspects ont été abordés : (i) la modélisation mathématique, à travers le développement du premier modèle tridimensionnel du réseau cérébral veineux et l'évaluation de l'importance des différentes hypothèses de modélisation par une approche de type analyse de sensibilité ; (ii) l'analyse théorique et numérique d'une nouvelle méthode de discrétisation, basée sur des multiplicateurs de Lagrange, qui a été développée afin de prendre en compte des conditions aux limites faisant intervenir la pression (méthode qui pourrait par ailleurs être intéressante dans le contexte plus général des réseaux hydrauliques) ; (iii) des contributions en calcul scientifique, avec l'objectif de développer un environnement de calcul haute performance, validé sur des cas-tests significatifs pour les applications en vue ; (iv) des avancées interdisciplinaires visant à créer des données angiographiques virtuelles à partir d'images médicales, dans le cadre du projet ANR VIVABRAIN. Une partie de ces résultats ont été obtenus au cours de la thèse de doctorat de Ranine Tarabay (Université de Strasbourg), co-encadrée avec Christophe Prud'homme (Université de Strasbourg) et Nicolas Passat (Université de Reims), soutenue en septembre 2016, ainsi que dans une collaboration avec Olivia Miraucourt (Université de Reims), encadrée par Stéphanie Salmon (Université de Reims), Hugues Talbot (ES-IEE Noisy Le Grand) et Nicolas Passat (Université de Reims), qui a soutenu sa thèse de doctorat en novembre 2016.

Le chapitre 2 aborde la question des modèles mathématiques et numériques réduits et multi-échelles pour les biofluides. Plus précisément, dans une première partie, nous avons développé plusieurs modèles 0d dans le but de comprendre la dynamique complexe des biofluides qui interagissent dans le système couplé oeil-cerveau. L'avantage de ces nouveaux modèles de type réseau est qu'ils fournissent une vue systémique, capable de capturer la dynamique globale de l'interaction entre le sang, le liquide céphalo-rachidien, les humeurs oculaires et le fluide interstitiel, dans l'oeil et dans le cerveau, tout en conservant une complexité mathématique raisonnable et un faible coût de calcul. Du point de vue clinique, les simulations ont montré des résultats intéressants dans le cadre de deux applications : l'étude de l'impact de la microgravité sur la vision des astronautes pendant les missions de longue durée dans l'espace et la comparaison de différentes stratégies thérapeutiques dans la prise en charge du glaucome. Dans la deuxième partie du chapitre, nous proposons un nouvel algorithme pour la résolution numérique d'un système couplé d'équations aux dérivées partielles et d'équations différentielles issu de la mécanique des fluides et motivé par des applications à la modélisation multi-échelles du flux sanguin à travers le système cardio-vasculaire. À nouveau, il convient de noter que le cadre conceptuel qui en résulte peut être pertinent dans le contexte plus général des réseaux hydrauliques. Les contributions présentées dans ce chapitre ont bénéficié des interactions au cours de la thèse de Lucia Carichino (IUPUI, soutenue en août 2016) et de celle de Simone Cassani (IUPUI, soutenue en août 2016), encadrées par Giovanna Guidoboni (Université du Missouri) et du travail accompli pendant la première année du doctorat de Lorenzo Sala (Université de Strasbourg), que je co-encadre actuellement avec Christophe Prud'homme (Université de Strasbourg) et Giovanna Guidoboni (Université du Missouri).

Dans le chapitre 3, la description de la dynamique des fluides donnée dans les chapitres précédents est enrichie pour prendre en compte les effets combinés de l'écoulement et de différentes structures, dans une perspective multi-physique. Une première contribution concerne la résolution numérique d'un problème couplé fluide-structure, décrivant la dynamique du flux sanguin par les équations de Navier-Stokes dans un domaine mobile, couplées avec les équations de l'élasticité linéaire qui régissent la déformation de la paroi du vaisseau. Ce modèle multi-physique couplé est ensuite utilisé pour simuler le flux sanguin dans l'aorte dans le but d'étudier le rôle de la géométrie de l'arche aortique dans les stades post-opératoires d'une malformation congénitale. Ces résultats font partie de la thèse de doctorat de Nicole Poussineau (UPMC, Paris 6), encadrée par Yvon Maday (UPMC, Paris 6) et soutenue en décembre 2007. Un autre type d'interaction fluide-structure, à savoir le cas d'une particule rigide évoluant passivement dans un fluide, est étudié dans la deuxième partie. L'accent est mis sur la nécessité de prendre en compte les contacts comme ingrédient clé de la simulation numérique directe des écoulements fluide-particules, notamment dans un régime dense. Nous avons proposé un nouveau schéma de discrétisation en temps qui fonctionne de manière robuste dans la situation difficile où la force de traînée (aussi connue comme la force de lubrification dans ce cas) tend vers l'infini très rapidement lorsque la particule s'approche du mur.

Mes recherches antérieures portaient sur l'étude de certaines questions qui apparaissent en théorie mathématique de l'élasticité, en utilisant des méthodes d'analyse mathématique et de géométrie différentielle. Au lieu de considérer la déformation du corps comme l'inconnue principale du

Introduction

This manuscript gathers my contributions focused on developing new mathematical and computational methods for analyzing biological flows as complex multiphysics and multiscale phenomena. The description of the underlying mechanisms stems from the basic principles of fluid dynamics and is translated into systems of partial or ordinary differential equations. The overall objective of this work is the study of these equations at the continuous and discrete levels, their coupling and the development of a reliable and efficient computational framework to implement various numerical methods to approximate the solutions to these problems. The numerical simulations incorporate realistic geometries, are thoroughly validated against experimental data and target specific biomedical applications.

The in silico approach, namely combining data-driven and physics-driven modeling in a unifying computational framework, has become increasingly popular in the domain of life sciences in the last years and has motivated the development of new mathematical and numerical concepts. Despite the significant advances in the in silico modeling of human physiology, understanding the complex dynamics governing the interplay between different fluids in the human body is still an extremely challenging field. Biological fluid description encompasses a wide range of spatial and temporal scales, from the molecular level to networks of a few meters, between a one-second heartbeat and a lifetime. The fluid dynamics is influenced by the interaction with surrounding tissues, that calls for a multiphysics approach. The geometrical representation is very complex and the availability of real data is scarce. Moreover, in experimental studies, multiple factors come into play and it is difficult to single out each contribution on the flow characteristics. In this direction, statistical analysis offers many techniques to help identify correlations among real data, even though limited insights can be gained on the cause-effect relationships beyond the correlations. In this perspective, mathematical and computational modeling approaches, based on the fundamental principles of physics and adapted to a biological framework, can help elucidate cause-effect relationships among interplaying factors and increase our understanding of these complex phenomena, in synergy with other approaches.

At this stage, an important point should be emphasized: what are the differences between modeling biological fluids and classical fluids? Compared to classical fluids, they serve crucial biological functions (e.g. transport of nutrients and oxygen, interaction with surrounding cells), that need to be accounted for in a mathematical model. In contrast with "inert" fluids, their properties change with age, disease, etc., as a consequence of biological processes. There are numerous challenges to study them in vivo, since most of the time biological fluids are difficult to access and study in their natural environment; in addition, they often show different properties and behaviors when examined in vitro or in vivo. Therefore, extending the classical fluid mechanics and mathematical theory to biological systems is not straightforward and represents an active field of research.

6

All these aspects result in (i) numerous challenging issues on how to set up a sound mathematical model, based on the biological and clinical literature; (ii) a compelling necessity to develop new theoretical approaches and numerical methods; (iii) an increasing need to design and implement flexible and efficient computational frameworks; (iv) important research efforts towards verification, validation and treatment of uncertainty. These aspects are crucial to guarantee that the model and its numerical solutions are meaningful from the biomedical viewpoint, thereby allowing the use of the resulting computational framework as a virtual laboratory where ideas can be tested and new hypotheses can be formulated.

In this context, the contributions described in this manuscript are a triple-point crossroads between mathematics, modeling and medicine. As summarized below and detailed in the next chapters, we will discuss the development of new mathematical models and numerical methods, illustrated by specific biomedical applications. Through the selected publications, each chapter presents the challenges and achievements on the theoretical issues and the practical aspects, and is concluded by a general discussion and perspectives. All these developments would not have been possible without a strong interaction with colleagues from mathematics, physics, computer science, biology and medicine. Their names appear in the list of references. I will particularly mention in this introduction the contributions of several PhD students who I had the opportunity to mentor, formally or informally, and develop joint publications.

The chapters are structured as follows:

In Chapter 1, we focus on three-dimensional models, in which the motion of a biofluid in a complex, realistic geometry is governed by the Navier-Stokes system in a domain with inlet and outlet boundaries. Different aspects were addressed: (i) mathematical modeling issues, through the development of the first three-dimensional model for the cerebral venous network and the assessment of the importance of different modeling assumptions through a sensitivity analysis approach; (ii) theoretical and numerical questions, concerning a new Lagrange multiplier-based numerical method accounting for boundary conditions involving pressure (that could be of interest in the more general context of hydraulic network-like systems); (iii) scientific computing contributions to the development of a high performance computing framework, validated on significant benchmarks and (iv) interdisciplinary investigations within the VIVABRAIN project, aiming at creating virtual angiographic data starting from real medical images. Part of these results were obtained during the PhD thesis of Ranine Tarabay (Univ. de Strasbourg), co-advised with Christophe Prud'homme (Univ. de Strasbourg) and Nicolas Passat (Univ. de Reims), who defended in September 2016 and in a collaboration involving Olivia Miraucourt, who defended her PhD thesis in November 2016 (Univ. de Reims), mentored by Stéphanie Salmon (Univ. de Reims), Hugues Talbot (ESIEE Noisy Le Grand) and Nicolas Passat (Univ. de Reims).

Chapter 2 deals with reduced and multiscale mathematical and computational models for biofluids. More precisely, in a first part we developed several 0d models with the aim of understanding the complex fluid dynamics in the coupled eye-cerebral system. The advantage of these new networkbased models is that they provide a systemic view, able to capture the overall dynamics of the interwoven physiology of blood, cerebrospinal fluid, ocular humors and interstitial fluids in the eye and in the brain, while maintaining a relatively accessible mathematical complexity and low computational costs. From the clinical viewpoint, simulations showed interesting features when applied to microgravity conditions and for therapeutical strategies in glaucoma management. In the second part of the chapter, we describe a new algorithm for the numerical solution of coupled systems of partial and ordinary differential equations for fluid flows, motivated by applications to blood flow modeling through the cardiovascular system from a multiscale perspective. Again, it should be noted that the resulting conceptual framework may be meaningful and applicable to a more general context of hydraulic networks. The contributions presented here benefit from the interactions during the PhD thesis of Lucia Carichino (IUPUI, defended in August 2016) and Simone Cassani (IUPUI, defended in August 2016), both mentored by Giovanna Guidoboni (Univ. of Missouri) and from the work achieved through the first year of the PhD of Lorenzo Sala (Univ. de Strasbourg), who I am currently co-advising with Christophe Prud'homme (Univ. de Strasbourg) and Giovanna Guidoboni (Univ. of Missouri).

In Chapter 3, the fluid dynamics description from the previous chapters is enriched to take into account the combined effects of flow and different structures, from a multiphysics perspective. A first contribution concerns the numerical solution of a coupled fluid-structure interaction model, based on a description of blood flow dynamics by the Navier-Stokes equations in a moving domain, coupled with the linear elasticity equations describing the vessel wall deformation. This coupled multiphysics model is subsequently utilized for simulating blood flow in the aorta with the aim of investigating the role of geometry of the aortic arch in the post-operative stages of a congenital malformation. These results are part of the PhD thesis of Nicole Poussineau (UPMC, Paris 6), defended in December 2007 and mentored by Yvon Maday (UPMC, Paris 6). Another type of fluidstructure interaction, namely the case of a rigid particle evolving passively in a fluid, is presented in the second part. The emphasis is the issue on how to handle contacts, as key ingredient in the direct numerical simulation of particulate flows, especially in the dense regime. We proposed a new timeadvancing scheme that works robustly in the challenging situation where the drag force (also known as the lubrication force in this case) tends to infinity very rapidly as the particle approaches the wall.

My earlier research dealt with the study of some questions that arise in the theory of elasticity, by using methods of mathematical analysis and differential geometry. Instead of considering the body deformation as the main unknown of the problem, the problem is expressed in terms of the geometrical properties of the elastic body under deformation. In this direction, we mention here references [START_REF] Szopos | On the recovery of a curve isometrically immersed in E n[END_REF] (and the abridged version [START_REF] Szopos | Sur la reconstruction d'une courbe isometriquement immergee dans un espace Euclidien[END_REF]), [START_REF] Szopos | On the recovery and continuity of a submanifold with boundary[END_REF] (and the abridged version [START_REF] Szopos | On the recovery and continuity of a submanifold with boundary in higher dimensions[END_REF]), and [12]. They contain results obtained during my PhD thesis and will not be developed in the sequel. For the sake of completeness, we also cite reference [START_REF] Szopos | An existence and uniqueness result for isometric immersions with little regularity[END_REF], providing a generalization of the results in [START_REF] Szopos | On the recovery and continuity of a submanifold with boundary[END_REF] obtained after my PhD thesis, but not described here.

Chapter 1

Three-dimensional models, analysis and efficient simulations for fluid equations. Applications to computational hemodynamics in the cerebral network

The goal of this chapter is to discuss the main challenges in blood flow mathematical and computational modeling in realistic geometries and to present our contributions within this framework. We focus here on the so-called macro-scale models, that are suitable when aiming for a threedimensional description of flow in large vessels, as derived in the case of large cerebral veins in Section 1.1. The fluid dynamics is described by the Navier-Stokes equations, endowed with suitable initial and boundary conditions. In Section 1.2, we emphasize the particular importance of boundary conditions and thoroughly analyze a novel formulation for the Stokes problem with non standard boundary conditions involving the pressure. The impact of different modeling assumptions on hemodynamic quantities of interest (velocity and wall shear stress) is investigated by means of a systematic numerical exploration in Section 1.3. From the computational viewpoint, large scale three-dimensional simulations are needed to account for the intrinsic 3d nature of these problems. We therefore developed a reliable and efficient computational framework, that we validate against experimental data, as described in Section 1.4. All these new contributions were integrated in an interdisciplinary research project aiming at creating virtual angiographic data starting from real medical images, and in the open-source software pipeline associated to this framework. Section 1.5 presents the challenges and achievements of this integration. Finally, further extensions and new directions inspired by the previous works are reviewed in Section 1.6.

Mathematical modeling of blood flow in the cerebral venous system

The circulatory system is a closed circuit through which blood flows and transports nutrients. Its main components are the arteries (blood vessels that carry blood away from the heart), the capillary bed (where the functional exchanges occur) and the veins (blood vessels that carry blood toward the heart). In the cerebral compartment, the vasculature has a very complex three-dimensional 10 CHAPTER 1. THREE DIMENSIONAL MODELS structure and blood flow is high, in order to properly satisfy brain's energy demands. Clinical data report values of global cerebral blood flow around 700 ml/min, accounting for more than 15% of total cardiac output [START_REF] Partington | Intracranial pressure and cerebral blood flow[END_REF], for only about 2% of the total body weight. Normal intracranial blood volume is around 150 ml, two thirds of which is in the venous system. Overall, these figures highlight the crucial role played by hemodynamical factors in the brain and the importance of understanding their influence in physiological and pathological situations.

A new three-dimensional model for large cerebral veins

The model derivation and the mathematical/fluid-mechanics study of blood flow in large cerebral arteries (e.g. carotid arteries, basilar trunk, and the circle of Willis) are the focus of multiple investigations, as reviewed for instance in [START_REF] Formaggia | of MS&A[END_REF] or very recently in [START_REF] Quarteroni | The cardiovascular system: Mathematical modelling, numerical algorithms and clinical applications[END_REF]. By contrast, the behavior of the venous cerebral system is rarely explored. This compartment has its own features, that are contributing to the difficulty of deriving an appropriate mathematical and computational model, among which:

• complex three-dimensional structure, see Figure 1.1, often asymmetric, presenting a considerably more variable pattern than the arterial system [START_REF] Schaller | Physiology of cerebral venous blood flow: from experimental data in animals to normal function in humans[END_REF];

• lack of parallelism between arterial and venous circulations, in contrast to the rest of the body [Sch04];

• highly individual variations of the venous outflow [START_REF] Stoquart-Elsankari | A phase-contrast MRI study of physiologic cerebral venous flow[END_REF].

Therefore, this subject has received far less attention than the arterial network and only few recent references are available in the literature to the best of our knowledge. A 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] Gadda | A new hemodynamic model for the study of cerebral venous outflow[END_REF] (lumped parameter model of the cerebral venous outflow), [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 [MT14a; MT14b] (closed-loop 1d model of the human circulation, with a special focus on the venous system). An alternative approach is to utilize three-dimensional (3d) models of a region of interest, in idealized or patient specific geometries. Classically, due to prohibitive computational costs, modeling is limited to relatively small samples of the vascular circuit, possibly coupled to 0d or 1d models. Except for the contributions presented here, we are aware of only one recent work [START_REF] Caiazzo | Computational haemodynamics in stenotic internal jugular veins[END_REF] proposing a 3d description of blood flow in rigid jugular veins coupled with a 1d model for major intracranial veins.

Our contribution [11] is to rely on 3d blood flow modeling in more complete vascular networks, in order to better describe blood flow patterns in complicated geometries. We proposed a complete computational modeling approach from medical image processing to numerical resolution. The realistic three-dimensional description of the geometry of the network, see Figure 1.1, 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 [10], see also Section 1.5. 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 model assumes that intracranial veins located between the skull and the brain are rigid.

Data used in the model are given in Table 1.1. However, since patient-specific blood flow information are usually not collected from routine clinical examinations, only scarce literature data are available for the velocity magnitude in cerebral veins, which adds new difficulties in deriving

ρ ∂u ∂t -2∇ • (µD(u)) + ρ(u • ∇)u + ∇p = 0, in Ω × I (1.1) ∇ • u = 0, in Ω × I (1.2)
endowed with appropriate initial and boundary conditions, where: u and p are the velocity and pressure of the fluid, D(u) = 1 2 (∇u + ∇u T ) is the strain rate tensor, σ(u, p) = -pI + 2µD(u) is the Cauchy stress tensor, ρ and µ are the density and dynamic viscosity of the fluid, respectively.

Concerning the initial status of the fluid velocity, it is well known that it has to be carefully prescribed, since it should be divergence-free to be admissible. Unfortunately, in hemodynamic computations, this quantity is usually unknown, hence chosen equal to zero everywhere or, as a better guess, as the solution of a stationary Stokes problem. We developed a solution to this problem in the preliminary work [START_REF] Salmon | Cerebral blood flow simulations in realistic geometries[END_REF].

The issue of boundary conditions is of primary importance in simulating blood flow and a huge literature has been dedicated to this topic in the last years, as reviewed for instance in [START_REF] Fouchet-Incaux | Artificial boundaries and formulations for the incompressible Navier-Stokes equations: applications to air and blood flows[END_REF][START_REF] Quarteroni | The cardiovascular system: Mathematical modelling, numerical algorithms and clinical applications[END_REF]. In the context of the cerebral venous network, we proposed an extensive discussion on this topic in [11], with the following partial conclusions:

• at the inflow, impose u = constant (small magnitude), since blood comes from the microcirculation, modeled by a quasi-steady/steady Stokes flow;

• at the wall, impose u = 0, since intracranial veins are constrained between a nearly incompressible brain and the rigid skull;

• at the outflow, impose in a first approximation σ(u, p)n = 0 (the "do-nothing" classical approach).

Perspectives. A lot of subsequent questions are triggered by these choices, as for instance: how small (in terms of order of magnitude) should the inflow velocity be? How a rigid vessel model (acceptable for intracranial veins) could be coupled with a deformable vessel model (that is more appropriate for jugular veins)? At the outflow: between using the "do-nothing" approach and coupling with a reduced model involving a lot of parameters that are unavailable, which choice is better? We postpone the answers to some of these questions to Section 1.2 for a different theoretical approach [START_REF] Bertoluzza | Boundary conditions involving pressure for the Stokes problem and applications in computational hemodynamics[END_REF] and to Section 1.3 for a sensitivity analysis study on the influence of different modeling assumptions [START_REF] Chabannes | Hemodynamic simulations in the cerebral venous network: A study on the influence of different modeling assumptions[END_REF]. We also discuss in Section 1.6 some questions left open and possible lines of approach.

Numerical results and discussion

The numerical strategy we implemented to compute approximate solutions to the previously derived model involves a time-scheme based on the characteristics method [START_REF] Pironneau | On the transport-diffusion algorithm and its applications to the Navier-Stokes equations[END_REF] and a spatial discretization of finite element type, with P 2 /P 1 inf-sup stable finite elements [START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF]. The numerical solving approach relies on the finite element library Freefem++ [START_REF] Hecht | New development in FreeFem++[END_REF]. More details about the parameters and numerical choices can be found in [11].

We display in Fig. 1.2 the simulated velocity field in the whole intracranial cerebral venous network. We were able to identify a complicated three-dimensional flow behavior, with noticeable recirculations in the confluence of the sinuses (bottom left panel). At the level of the superior sagittal sinus, the flow is laminar (bottom right panel). An asymmetric behavior appears in a recirculation zone in the entrance segment of the right transverse sinus. This could be explained, at least partly, by the asymmetric architecture of the venous network. An important question is the validation of these results, that we investigated as follows: 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 [Sch04] and 15.2 ± 3 • 10 -2 m/s [Gid+96], 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 phase-contrast 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 [START_REF] Müller | A global multiscale mathematical model for the human circulation with emphasis on the venous system[END_REF][START_REF] Müller | Enhanced global mathematical model for studying cerebral venous blood flow[END_REF]. Comparisons between the model predictions of flow rate values at selected cross-sections with MRI measurements reported in the literature are presented in Fig. 1.3.

Although these comparisons show satisfactory results, it should be noted that the clinically indicated patient measurements are not taken on the same geometry as the one used to perform our simulations, but retrieved from reported data in the literature. Therefore, improving the validation process, using either in vitro or in vivo data is still a challenge, that will be further explored in Section 1.4.

A novel formulation of the Stokes system involving pressure boundary conditions

The issue of boundary conditions when modeling blood flow in the circulatory system is of major importance and matter of intense research. From the mathematical viewpoint, the well-posedness Section 3], see also [START_REF] Fouchet-Incaux | Artificial boundaries and formulations for the incompressible Navier-Stokes equations: applications to air and blood flows[END_REF]. When performing three-dimensional simulations for blood flow modeling, the domain is reduced to a region of interest and therefore, its boundary is composed of two parts: a physical boundary (corresponding to the vessel wall) and an artificial boundary (at the level where the vessel is truncated). On one hand, at the vessel wall, the no-slip Dirichlet condition is natural for a viscous fluid contained in a rigid domain, since the viscous effect constraints the fluid particles to adhere to the wall. On the other hand, at the artificial boundaries, different formulations with boundary conditions involving components of the velocity field, stresses or pressure are of interest. Indeed, they should be able to take into account the rest of the closed circuit representing the circulatory system. Moreover, in order to achieve physiological simulations, they should derive from clinical measurements. However, in this case they might give rise to so-called defective boundary conditions: usually only average data are available to "feed" mathematical and computational models that require pointwise data instead, see [START_REF] Formaggia | Prescription of general defective boundary conditions in fluid-dynamics[END_REF] or [START_REF] Porpora | Numerical treatment of boundary conditions to replace lateral branches in hemodynamics[END_REF] for different strategies in the context of blood flow modeling.

In the particular case of the cerebral venous network, one of the main mechanisms that drives the flow and should be taken into account through appropriate boundary conditions is the pressure drop. More generally, when modeling hydraulic network-like systems, for instance oil ducts, water supply, microfluidic channels or biological flows, we are interested in the case when the velocity field is imposed on one part of the boundary and pressure values are prescribed, together with the condition of no tangential flow, on the remaining part. We present in this section a new method we proposed in [START_REF] Bertoluzza | Boundary conditions involving pressure for the Stokes problem and applications in computational hemodynamics[END_REF] to take into account these non-standard boundary conditions, both at the continuous and the discrete levels, in a finite-element framework.

General framework and the Lagrange multiplier formulation

We consider the steady state of a viscous incompressible fluid at low Reynolds number, described by the velocity and pressure fields u and p that satisfy the following Stokes equations:

-2µ∇ • (D(u)) + ∇p = ρf , in Ω, (1.3) ∇ • u = 0, in Ω, (1.4)
u = 0, on Γ 1 , (1.5)
u × n = 0, on Γ 2 , and (1.6)
p = p 0 , on Γ 2 , (1.7)
where A variational formulation taking into account this type of boundary conditions was first introduced in the seminal works [START_REF] Pironneau | Conditions aux limites sur la pression pour les équations de Stokes et de Navier-Stokes[END_REF][START_REF] Bègue | A nouveau sur les équations de Stokes et de Navier-Stokes avec des conditions aux limites sur la pression[END_REF][START_REF] Conca | The Stokes and Navier-Stokes equations with boundary conditions involving the pressure[END_REF]. A lot of subsequent literature was devoted to this topic, see [START_REF] Bertoluzza | Boundary conditions involving pressure for the Stokes problem and applications in computational hemodynamics[END_REF] for an in depth discussion about the method in the context of the existing literature.

∂Ω = Γ1 ∪ Γ2 , with Γ 1 ∩ Γ 2 = ∅
Remark 1.2.1. Previous works [Pir86; Bèg+87; CMP94; Con+95; BCRY15] classically take into account non standard boundary conditions of type (1.5-1.7) by expressing the conservation of the momentum in terms of the Laplacian of the velocity and then using as a key ingredient the rotational formulation for the equation, based on:

-Δu = ∇ × (∇ × u) -∇(∇ • u).
Although at a continuous level the two formulations are equivalent, since

∇ • u = 0 ⇒ ∇ • (∇u + ∇u T ) = Δu,
from a modeling standpoint it may be useful to work with the symmetric tensor. For instance, in fluid-structure problems, formulation (1.3-1.7) gives directly the natural boundary condition for the structure problem in terms of the force exerted by the fluid on its boundary. We thus focus hereafter on the formulation in terms of the divergence of the symmetric gradient (1.3).

CHAPTER 1. THREE DIMENSIONAL MODELS

We showed in [START_REF] Bertoluzza | Boundary conditions involving pressure for the Stokes problem and applications in computational hemodynamics[END_REF] that the problem can be written as follows:

Problem 1.2.2. Find u ∈ V , p ∈ M , λ ∈ Λ such that for all v ∈ V , q ∈ M , η ∈ Λ 2µ � Ω D(u) : D(v) dx - � Ω p∇ • v dx -c(v, λ) = ρ � Ω f • v dx - � Γ 2 p 0 n • v ds (1.8) � Ω q∇ • u dx = 0 (1.9) c(u, η) = 0. (1.10)
The functional spaces are defined as

V = {v ∈ [H 1 (Ω)] d : v = 0 on Γ 1 }, M = L 2 (Ω), and Λ = [H -1/2 (Γ 2 )] d-1 . The bilinear form c : V × Λ → R is given by c(u, λ) = � Γ 2 u • i(λ) dx where i : [H -1/2 (Γ 2 )] d-1 → T is an isomorphism and T = � ζ ∈ [H -1/2 (Γ 2 )] d : ζ • n = 0 � .
Details on how the bilinear form c is built in practice in two and three dimensions, by making explicit the isomorphism between T and Λ can be found in [START_REF] Bertoluzza | Boundary conditions involving pressure for the Stokes problem and applications in computational hemodynamics[END_REF]. We have the following theorem:

Theorem 1.2.3. Problem 1.2.2 admits a unique solution (u, p, λ) which verifies Discretization of Problem 1.2.2. We introduce a compatible tessellation T h of the domain Ω in tetrahedral or hexahedral elements [START_REF] Bertoluzza | Boundary conditions involving pressure for the Stokes problem and applications in computational hemodynamics[END_REF]. On T h , we introduce piecewise polynomial spaces V h ⊆ V , Q h ⊂ M , respectively approximating velocity and pressure, and we assume that such spaces satisfy the standard inf-sup condition inf

�u� 1,Ω + �p� 0,Ω + �λ� -1/2,Γ 2 � �f � V � + inf v∈V � Γ 2 p 0 n • v ds �v� 1,Ω � �f � 0,Ω + �p 0 � 0,Γ 2 . Moreover, if u ∈ [C 2 (Ω)] d , p ∈ C 1 (Ω),
p h ∈Q 0 h sup u h ∈V h ∩[H 1 0 (Ω)] 3 � Ω p h ∇ • u h dx �u h � 1,Ω �p h � 0,Ω � 1. (1.11) (where Q 0 h = � q h ∈ Q h : � Ω q h = 0 �
), so that they provide a stable discretization of the Stokes problem with standard Dirichlet boundary conditions. We now observe that 3 where W h is itself a finite element space on the two-dimensional mesh T Γ 2 h induced on Γ 2 by the threedimensional tessellation T h . Remark (see the definition of the space V ) that the functions in W h satisfy homogeneous boundary conditions on ∂Γ 2 . We then let Λ h = [W h ] 2 . We consider the following discrete problem:

V h | Γ 2 = [W h ]
Problem 1.2.5. Find u h ∈ V h , p h ∈ Q h , λ h ∈ Λ h such that for all v h ∈ V h , q h ∈ Q h , η h ∈ Λ h 2µ � Ω D(u h ) : D(v h ) dx - � Ω p h ∇ • v h dx -c(v h , λ h ) = ρ � Ω f • v h dx - � Γ 2 p 0 n • v h ds, (1.12) � Ω q h ∇ • u h dx = 0, (1.13) c(u h , η h ) = 0. (1.14)
The following theorem then holds:

Theorem 1.2.6. There exists h 0 such that, if h ≤ h 0 , Problem 1.2.5 admits a unique solution (u h , p h , λ h ) which verifies

�u h � 1,Ω + �p h � 0,Ω + �λ h � -1/2,Γ 2 � �f � 0,Ω + �p 0 � 0,Γ 2 .
Moreover the following error estimate holds:

�u -u h � 1,Ω + �p -p h � 0,Ω � inf v h ∈V 0 h �u -v h � 1,Ω + inf q h ∈Q h �q -q h � 0,Ω , where V 0 h = {u h ∈ V h : c(u h , λ h ) = 0, ∀λ h ∈ Λ h } .
Once again, as in the continuous case, the proof of Theorem 1.2.6 reduces to prove two inf-sup conditions. A direct consequence of this result is that for the classical case of Taylor-Hood inf-sup stable finite element spaces [START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF]:

V h = {u ∈ [C 0 (Ω)] 3 : ∀K ∈ T h u| K ∈ [P k (K)] 3 }, (1.15)
Q h = {p ∈ C 0 (Ω) : ∀K ∈ T h p| K ∈ P k-1 (K)}. (1.16) and for u ∈ [H k+1 (Ω)] 3 and p ∈ H k (Ω) we have �u -u h � 1,Ω + �p -p h � 0,Ω � h k (�u� k+1,Ω + �p� k,Ω ). (1.17)
We thus expect optimal convergence rates provided the solution has sufficient regularity.

Remark 1.2.7. We emphasize that there is no reason why the multiplier i(λ) = τ (u, p) should vanish at the boundary of Γ 2 . Therefore, the proposed discretization cannot, in general, yield an optimal approximation of the Lagrange multiplier. Nevertheless, since V 0 h ⊂ V 0 , the approximation properties for the Lagrange multiplier do not enter the error estimate in Theorem 1.2.6, and we get an optimal error estimate for both velocity and pressure.

Remark 1.2.8. Throughout the study, we assumed that Γ 2 was a flat surface (or, more precisely, we assumed that n was constant on the connected components of Γ 2 ). Let us consider two cases in which this assumption is not satisfied.

If a connected component of Γ 2 is the union of two (or more) flat subregions sharing a vertex (in 2d) or an edge (in 3d), we observe that, at the continuous level, if n 1 and n 2 are constant unit normals to the two subregions with direction chosen in such a way that on the common vertex or edge we have u = |u|n 1 = |u|n 2 , if n 1 � = n 2 then u = 0, so that the solution vanishes at the interface between the two subregions. At the discrete level one needs then to strongly force the velocity to vanish on such interface. Once this is done, the analysis presented above remains valid.

If, on the other hand, Γ 2 is a curved surface, the situation is more complex. At the continuous level, we show that the natural boundary condition implicit in equation (1.8) is not (1.7), but rather

p + 2µ|u|κ = p 0 , on Γ 2 , (1.18)
where κ is the mean curvature of Γ 2 . Nevertheless, Problem 1.2.2 is still well posed, and Theorem 1.2.3 still holds, provided equation (1.7) is replaced by equation (1.18). Things are more complex when it comes to discretization Problem 1.2.2, since part of the arguments used in the proof of Theorem 1.2.6 do not hold for curved boundaries. In addition, if the normal to the discrete boundary has jumps (which would automatically happen when approximating a curved boundary with a finite element mesh), similar arguments show that the whole method would then be non conforming. Remark that we might also need to resort to a non conforming discretization if we drop the requirement that the tessellation T h is compatible with the splitting of ∂Ω into Γ 1 ∪ Γ 2 .

Numerical results and applications

The computational framework is developed in Feel++, Finite Element Embedded Library in C++ [Fee], see for more details Section 1.4. Regarding the specific implementation of the proposed methodology, we would like to point out some non-standard aspects, namely the treatment of the terms associated to the Lagrange multipliers. Feel++ provides a mechanism to extract submeshes of faces and keep a relation between the extracted mesh and the parent mesh. The relation is necessary to ensure an efficient treatment of the coupling terms between the velocity and the Lagrange multipliers. The geometrical data, i.e the normals, are automatically deduced from the parent mesh.

We evaluated the performances of the method proposed in Section 1.2.1 when solving two types of problems: (i) study of the convergence properties of the method for different choices of finite elements, (results not reported here, available in [START_REF] Bertoluzza | Boundary conditions involving pressure for the Stokes problem and applications in computational hemodynamics[END_REF]); (ii) a 3d computational model of the cerebral venous blood flow similar to the one described in Section 1.1, but incorporating pressure boundary conditions (in a first approximation for low Reynolds numbers).

Figure 1.4 displays the pressure field (top left panel) and instantaneous streamlines, colored with velocity magnitude (top right panel), illustrating the pressure drop effect and a complicated threedimensional flow behavior. The overall dynamics shares common features with the one we found in Section 1.1, but a direct comparison would not be meaningful at this stage, since here only the Stokes equations are solved. The order of magnitude of the maximum velocity is slightly higher than values retrieved in the clinical literature, see for instance [START_REF] Schaller | Physiology of cerebral venous blood flow: from experimental data in animals to normal function in humans[END_REF], therefore more precise values need to be included in further work. However, the development of a computational model able to capture, to this level of accuracy, different features of the flow can be seen as a very promising approach for analyzing, by means of numerical simulations, the dynamics of flow patterns in morphologically complex vascular districts.

Remark 1.2.9. A zoom on some inlet, respectively outlet sections is presented in Figure 1.5, demonstrating that the flow is normal to both inflow and outflow surfaces. We highlight the interest of imposing the pressure value and the zero tangential component of the velocity in this context: the current formulation allows to retrieve a Poiseuille-like behavior that is physically meaningful when dealing with artificial boundary conditions, while keeping the viscous stress tensor in the expression of the Stokes problem, useful from a modeling standpoint. In contrast, as noted in [LG81; HRT96], using the symmetric gradient 1 2 (∇u + ∇u T ) and prescribing the normal stress at the outlet lead to a non-physical representation of the flow: the velocity vectors "spread" like at the end of a pipe, instead of mimicking the fact that the network continues after this artificial section. Alternatively, the non symmetric tensor ∇u can be used to recover the Poiseuille exact solution in a cylinder, but the physical meaning of such a boundary condition is not clear. Perspectives. The current methodology should be further developed, in particular by (i) devising an adapted discretization strategy for the case of a curved boundary, in order to overcome the difficulties briefly discussed in Section 1.

Influence of different modeling assumptions

The complexity of the human circulatory tree and the underlying bio-mechanical phenomena call for various simplifications in view of designing and implementing mathematical and computational models. The relative importance of these assumptions is difficult to establish and quantify. Therefore a lot of open questions still remain, see Figure 1.6 and extended discussion in [START_REF] Steinman | Assumptions in modelling of large artery hemodynamics[END_REF]. In addition, as already emphasized in Section 1.1, the venous network is far less known and studied in comparison to the arterial one. For arteries, recent contributions assessed the impact of viscosity models and flow conditions in aneurysms [START_REF] Evju | A study of wall shear stress in 12 aneurysms with respect to different viscosity models and flow conditions[END_REF] or of the assumption of laminar flow in computational hemodynamics, taking cerebral aneurysms as an illustrative example [START_REF] Evju | On the assumption of laminar flow in physiological flows: Cerebral aneurysms as an illustrative example[END_REF].

In line with these studies, the goal of [START_REF] Chabannes | Hemodynamic simulations in the cerebral venous network: A study on the influence of different modeling assumptions[END_REF], summarized in this section, is to assess in a sound mathematical and computational framework the effect of different modeling assumptions on cerebral venous blood flow dynamics at a macroscopic scale.

Sensitivity analysis framework

Formulation of the problem. We consider, as a starting point, the Navier-Stokes system (1.1)-(1.2) for large and medium-sized cerebral veins introduced in Section 1.1.1 with two extensions:

• a Generalized Newtonian constitutive law for blood flow, in addition to the Newtonian model from Section 1.1.1;

• the coupling with a simplified 0d model at the outflow, as an alternative to the traction-free boundary condition used in Section 1.1.1.

More precisely, the viscosity is considered either constant, when adopting a Newtonian constitutive model, or as a function of the shear rate:

γ = � 2tr � D (u) 2 � , (1.19)
when using a Generalized Newtonian model. In [START_REF] Chabannes | Hemodynamic simulations in the cerebral venous network: A study on the influence of different modeling assumptions[END_REF], we focused on the comparison between the Newtonian model, and the Carreau and Carreau-Yasuda Generalized Newtonian constitutive models for blood [START_REF] Formaggia | of MS&A[END_REF]Chap. 6], that return the viscosity as a function of the shear rate by the following equation:

µ( γ) -µ ∞ µ 0 -µ ∞ = (1 + (λ γ) a ) n-1 2 , (1.20)
where λ is a time constant, a and n are dimensionless parameters used to differentiate between the two models, and µ 0 and µ ∞ are the viscosities at zero and infinite shear rate, respectively. Regarding the boundary conditions: at the inflow, we impose a constant profile of small magnitude and at the vessel wall u = 0, since we considered it to be rigid. A more general setting is used for the outflow: boundary conditions are prescribed either by using a traction-free condition as in Section 1.1.1:

σ(u, p)n = 0 on Γ out , (1.21)
or by introducing the coupling with a three-element Windkessel model [FQV09, Chap. 10], in order to take into account the downstream vasculature. In the latter case, the condition reads:

σ(u, p)n = -P l n on Γ out , (1.22)
where P l (the proximal pressure) is obtained by solving

   C d,l dπ l dt + π l R d,l = Q l P l = R p,l Q l + π l , (1.23)
for a given value of the flux on the outlet

Q l = � Γout u • n dx.
Discretization. We use a fully implicit time discretization using BDF η scheme, including for the non-Newtonian models, with non-linear solving handled by a Newton method. The timediscretization of Equations (1.1)-(1.2) in the Newtonian case is written as follows:

ρ � � η k=0 α k u n+1-k Δt + u n+1 • ∇u n+1 � -∇ • (-p n+1 I + 2µD(u n+1 )) = 0, ∇ • (u n+1 ) = 0, (1.24)
where (α k ) k=0,η are the BDF η coefficients for the time derivative of u and the subscript η refers to the order of the scheme. We adopted the BDF 1 (η = 1) and BDF 2 (η = 2) schemes for the time approximation. The spatial discretization is handled via a inf-sup stable finite element (Taylor-Hood) P 2 /P 1 [START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF]. Our numerical strategy is in line with [START_REF] Valen-Sendstad | Mind the gap: impact of computational fluid dynamics solution strategy on prediction of intracranial aneurysm hemodynamics and rupture status indicators[END_REF] as we performed high resolution simulations and not normal resolution simulations: both the spatial and the temporal discretization methods are second-order schemes and a small time step (Δt = 10 -3 s) was chosen, in order to adequately resolve the complex flow features.

Sensitivity analysis framework

The idea stems from [START_REF] Evju | Sensitivity analysis of simulated blood flow in cerebral aneurysms[END_REF], where the author defines two metrics, one for the velocity and one for the wall shear stress, in order to improve the visual comparison between snapshots of solutions by assessing the differences between two solutions in a more quantitative manner. Both metrics require to define a computed reference solution and then the computation of the time evolution of the average on the entire mesh of some specific quantities, which reflect point-wise differences between velocities and the wall shear stresses, respectively. This

CHAPTER 1. THREE DIMENSIONAL MODELS

approach has in addition the advantage of assessing the time evolution of these differences, whereas comparisons of snapshots are only showing spatial differences at some selected instants.

The first metric measures the difference between the computed solution (u C ) and the computed reference solution (u R ). To this end, following [START_REF] Evju | Sensitivity analysis of simulated blood flow in cerebral aneurysms[END_REF], we define the spatial metric for each time step

S u R (u C ) = 1 -βe -a -γe -m , (1.25)
where

a = 1 π cos -1 � u R • u C �u C ��u C � � , 0 ≤ a ≤ 1, (1.26) m = �u R -u C � �u R � , 0 ≤ m and β = γ = 1 2 . (1.27)
The value a corresponds to the scaled angle between the two vectors and the value m to relative difference in their order of magnitude. The spatial average is then computed by:

D u = 1 V � Ω S u R (u C ) dx, where V = � Ω 1 dx. (1.28)
The second metric measures the difference between the computed reference solution of the wall shear stress τ R and another numerical solution of the wall shear stress τ C . Again, following [START_REF] Evju | Sensitivity analysis of simulated blood flow in cerebral aneurysms[END_REF], we define the spatial metric for each time step as

T τ R (τ C ) = 1 -e -t , t = � � � � τ R -τ C τ R � � � � (1.29)
and its spatial average

D τ = 1 S � ∂Ω T τ R (τ C ) dx, where S = � ∂Ω 1 dx. (1.30) Remark 1.3.1.
In practice, very intensive parallel computation are required to obtain the metrics (the studies used from 64 to 512 cores and ranged from 1 million to about 10 millions degrees of freedom). In our case, they were almost as expensive as some of the actual numerical simulations and therefore an efficient high performance computational framework was needed, see for more details Section 1.4.

Numerical results and discussion

The sensitivity analysis study was carried out in the open-source library Feel++ [Fee]. First, we carefully assessed the impact of the numerical strategy in use, see [START_REF] Chabannes | Hemodynamic simulations in the cerebral venous network: A study on the influence of different modeling assumptions[END_REF] and second, we investigated the impact of several modeling assumptions in a controlled numerical environment.

We illustrate here the outcomes of our study in the specific case of the influence of the inlet velocity magnitude imposed at each inflow. Detailed results for the other modeling assumptions (outflow treatment and viscosity models) are reported in [START_REF] Chabannes | Hemodynamic simulations in the cerebral venous network: A study on the influence of different modeling assumptions[END_REF]. The simulations were performed for several models (Newtonian 3d, non-Newtonian 3d and Newtonian 3d-0d). In each case, we take the reference solution with an inlet velocity magnitude equal to 10 -3 m s . Then, we compute the influence of each model by taking v in = 20 -3 m s and v in = 30 -3 m s , respectively. The results are displayed in Figure 1.7, where a grid refinement study is reported, and in Figure 1.8, where different models are compared. The mesh convergence tests were performed on three meshes M 0, M 1 and M 2 described in Table 1.2 and applied to the computed reference solution u R . Both metrics are persistently higher in comparison with the case of the outflow or rheological model, with a significant non-linear increase in the case of the wall shear stress metric. This should be compared to a Poiseuille approximation, frequently used in blood flow modeling, that predicts a linear increase of the wall shear stress when the velocity increases linearly. In the present case, where blood flow dynamics is complex and the geometry is also very complicated, this approximation is no more valid. It can also be noted that Newtonian 3d and Newtonian 3d-0d models present the same sensitivity for both velocity and wall shear stress metric; the same behavior can be found when comparing the two non-Newtonian models. Moreover, we can see in these figures that the non-Newtonian models have less influence compared to the Newtonian models on the same quantities. This effect is more significant in wall shear stress metric, but it is also visible in the velocity metric. Conclusions and outlook. The results of the present study showed that for cerebral veins blood flow modeling, the impact of setting the inlet boundary condition on the forces created by blood flow, is likely greater than for other modeling assumptions. Therefore, they highlighted the importance of deriving values for these conditions from clinically measured data at some probe locations, in order to enhance the accuracy of the computed hemodynamical quantities of interest. These findings are significant in the perspective of the integration of the computational modeling step in a full pipeline, where the interaction at the stage of acquisition of data and also in the validation process takes place, see Section 1.5. From the mathematical and numerical standpoint, only a deterministic approach was used in the present sensitivity analysis study, but next steps should also include quantification of statistical variability in these data. In this direction, we proposed a first attempt for a simplified 0d model in [START_REF] Szopos | Mathematical modeling of aqueous humor flow and intraocular pressure under uncertainty: towards individualized glaucoma management[END_REF], see Section 2.2.1.

Large-scale blood flow simulations and validation

In this section we present our continuous efforts to develop a powerful and flexible computational framework [START_REF] Caldini-Queiros | Towards large-scale three-dimensional blood flow simulations in realistic geometries[END_REF] for numerically solving problems arising in hemodynamics, such as [4; 1]. We also discuss our contribution towards improving the reliability and reproducibility of computational studies by performing a thorough validation of the fluid solver against experimental data [START_REF] Chabannes | High order finite element simulations for fluid dynamics validated by experimental data from the FDA benchmark nozzle model[END_REF].

Computational framework

The core of the framework is the Finite Element Embedded Library in C++ Feel++, that allows to use a very wide range of Galerkin methods and advanced numerical techniques such as domain decomposition (including mortar and three fields methods), fictitious domain or certified reduced basis. The ingredients include a very expressive embedded language, seamless interpolation, mesh adaption and seamless parallelization. Feel++ provides a mathematical kernel for solving partial differential equation using arbitrary order Galerkin methods in 1d, 2d, 3d and on manifolds using simplices and hypercubes meshes [Pru+12]:

• a polynomial library allowing for a wide range polynomial expansions including H div and H curl elements;

• a light interface to Boost.UBlas, Eigen3 and PETSc [Bal+16]/SLEPc as well as a scalable in-house solution strategy;

• a language for Galerkin methods starting with fundamental concepts such as function spaces, forms, operators, functionals and integrals;

• a framework that allows user codes to scale seamlessly from single core computation to thousands of cores and enables hybrid computing.

First steps towards large-scale three-dimensional blood flow simulations in realistic geometries were achieved in [START_REF] Caldini-Queiros | Towards large-scale three-dimensional blood flow simulations in realistic geometries[END_REF], where we investigated several issues: (i) handling various boundary conditions settings allowing for a flexible framework with respect to the type of input data (velocity, pressure, flow rate . . .); (ii) handling of the discretization errors not only with respect to the physical fields (velocity and pressure) but also with respect to the geometry; (iii) dealing with the associated large computational cost, requiring high performance computing, through strong and weak scalability studies.

On the basis of the problems studied in [4; 1], we added new capabilities in the algebraic solving framework. To give an insight about the importance of the preconditioning strategy when solving complex flow problems, we gather in Table 1.3 results allowing for a direct comparison between three possible choices in term of preconditioners when computing the solution of a pressure-driven Stokes flow in the cerebral venous network [START_REF] Bertoluzza | Boundary conditions involving pressure for the Stokes problem and applications in computational hemodynamics[END_REF] , the number of iterations grows strongly with respect to the problem size. Further refinements regarding the different choices are required and will be subject of future research.

The FDA benchmark nozzle model

A challenging benchmark was proposed by the US Food and Drug Administration (FDA) in [START_REF] Hariharan | Multilaboratory particle image velocimetry analysis of the FDA benchmark nozzle model to support validation of computational fluid dynamics simulations[END_REF] in order to assess the stability, accuracy and robustness of computational methods in different physiological regimes. The findings of 28 blinded investigations were reported in [START_REF] Stewart | Results of FDAs first interlaboratory computational study of a nozzle with a sudden contraction and conical diffuser[END_REF] and, as critically analyzed in [START_REF] Sotiropoulos | Computational fluid dynamics for medical device design and evaluation: are we there yet?[END_REF], practically all CFD solvers failed to predict results that agreed in a satisfactory manner with the experimental data. Several subsequent papers tackled this question, by employing different numerical approaches, see [START_REF] Chabannes | High order finite element simulations for fluid dynamics validated by experimental data from the FDA benchmark nozzle model[END_REF] for a more detailed discussion.

The benchmark provides a comprehensive dataset of experimental measures using a well-defined geometry corresponding to an idealized medical device (see Figure 1.9 for a schematic sketch of the domain and [Har+11, Sec. 2.1] for more details). Five sets of data spanning laminar, transitional and turbulent regimes are made available on-line and the device was designed to feature accelerating, decelerating and recirculating flow (all of which occur in real medical devices). We investigated in CHAPTER 1. THREE DIMENSIONAL MODELS Figure 1.9: FDA nozzle sketch and specifications (from [START_REF] Hariharan | Multilaboratory particle image velocimetry analysis of the FDA benchmark nozzle model to support validation of computational fluid dynamics simulations[END_REF]).

[5] three flow regimes, corresponding to a Reynold number at the throat level Re T hroat = 500, 2000 and 3500, by using a direct numerical simulation method for the Navier-Stokes equations. A detailed report of the numerical method and parameters used for each test can be found in [START_REF] Chabannes | High order finite element simulations for fluid dynamics validated by experimental data from the FDA benchmark nozzle model[END_REF]. In particular we implemented and compared low order as well as high order approximations including for the geometry and we discussed some issues not previously reported in the literature.

The comparison with experimental data is made in terms of (i) wall pressure difference (normalized to the mean throat velocity) versus axial distance; and (ii) axial component of the velocity (normalized to the mean inlet velocity) along the centerline:

Δp norm = p z -p z=0 1 2 ρ f u 2 t and u norm z = u z u i , where u t = 4Q πD 2 t , u i = 4Q πD 2 i , (1.31)
and Q is the volumetric flow rate. Furthermore, two validation metrics reported in [START_REF] Stewart | Results of FDAs first interlaboratory computational study of a nozzle with a sudden contraction and conical diffuser[END_REF] were computed: a conservation of mass error metric E Q (on a percentage basis) and a general validation metric E z comparing average experimental velocity data with computed axial velocities.

We only provide here illustrative numerical results corresponding to Re T hroat = 500; for all the other cases, see [START_REF] Chabannes | High order finite element simulations for fluid dynamics validated by experimental data from the FDA benchmark nozzle model[END_REF]. Four mesh refinement levels, denoted M0-M3, from the coarsest to the finest level, were constructed. Several polynomial order approximations were used and the notation P N +1 P N G kgeo specifies the discretization spaces for the velocity, pressure, and geometry, respectively. Figure 1.10 shows the results for the normalized axial velocity and the normalized pressure difference along the z axis, respectively. In each case, we can see very satisfactory agreement with the experimental data.

We illustrate in Figure 1.11 the computation of metrics E z and E Q for several mesh refinements at Re T hroat = 500. The metric E z takes small values in each numerical experiment, identifying a good agreement between computed and experimental data, and displays only small variations with respect to mesh refinement. On the other hand, the metric E Q is more sensitive to this factor: error does not exceed the ∼ 2%, except for the coarse mesh M 0 where, in two locations, the error increases up to ∼ 10%. Interestingly, we note that the P 3 P 2 G 1 approximation does not improve the results for the coarse mesh, but that a satisfactory error below 2% is retrieved when using a

P 2 P 1 G 2 approximation.
Additional tests to complement the study of the impact of high order approximation are ongoing and an extended version of [START_REF] Chabannes | High order finite element simulations for fluid dynamics validated by experimental data from the FDA benchmark nozzle model[END_REF] is planned to be submitted for publication.

An open-source framework to generate virtual MRA images from real MRA images

In the previous sections, we presented several contributions to develop a sound mathematical and computational framework for hemodynamic simulations, with an emphasis on applications to the cerebral venous network. This methodology was integrated in the interdisciplinary program VIVABRAIN Virtual angiography simulation from 3D and 3D+t brain vascular models [Viv], that aims at starting from real MR angiographic data to finally lead to the generation of virtual MR angiographic data. We summarize hereafter the main steps and challenges of these developments.

As explained in Section 1.1, the realistic data we used as a starting point were acquired using different Magnetic Resonance Angiography (MRA) sequences. MRA is an image acquisition technique particularly used, for instance, in the case of detection and diagnosis of vascular alterations (e.g., stenoses, aneurysms, thromboses), and patient follow-up after treatment (e.g., by stents, coils). MRA is also an effective tool for clinical research devoted to the understanding of human macro and mesovascular physiology, since it provides non-invasive/non-ionizing modalities (e.g., Phase Contrast (PC) or Time-of-Flight (TOF) MRA), available in 2d, 3d, and even 3d+time.

Due to the nature of structures visualized in MRA, namely thin, elongated, curvy vessels, the development of specific image analysis tools has been an active research area during the last decades [START_REF] Lesage | A review of 3D vessel lumen segmentation techniques: Models, features and extraction schemes[END_REF]. However, contrary to other morphological structures (e.g., brain tissues and structures, with BrainWeb [START_REF] Cocosco | Brainweb: Online interface to a 3D MRI simulated brain database[END_REF]), there is no common framework to facilitate the development and validation of related image analysis methods. In particular, despite a few attempts (such as [START_REF] Schaap | Standardized evaluation methodology and reference database for evaluating coronary artery centerline extraction algorithms[END_REF]), that rely on datasets equipped with manual ground truth, there is no vascular analogue of Brain-Web, i.e., a framework for generating efficiently virtual MRA images naturally equipped with a ground-truth and/or associated to a real MRA. Stemming on this fact, we have been working, in the context of the VIVABRAIN project with the final goal to tackle this issue, knowing that the case of virtual MRA was not intensively considered, in particular via the efficient coupling of CFD and MRI simulation.

The framework relies on a pipeline composed of five main steps, going from real MRA to virtual MRA, see Figure 1.12:

(i) MRA acquisition;

(ii) segmentation of vascular volumes from such data and design of a vascular model, in particular by explicitly modeling anatomical and physiological features;

(iii) construction of computational meshes of the vascular structures;

(iv) simulation of the flowing blood in the designed vascular model, leading to velocity and pressure fields;

(v) generation of virtual MRA images based on the previously computed information.

Our contributions are inserted in step (iv), with natural ramifications in the pre-and post-processing steps, namely building the computational meshes before the simulations and treating the results afterwards. Besides the inherent challenges related to the mathematical and computational approach, additional difficulties come into play. In particular, the complexity of the data to handle at the interface between the steps and their representation is a crucial point: for instance, fluid simulations are performed in an Eulerian framework, whereas MRA simulation in performed by These research efforts were first presented in [10], where the workflow from medical images to numerical simulations was discussed. A significant step forward was achieved by junior researchers of the VIVABRAIN project in [START_REF] Ancel | PHANTOM project: development and validation of the pipeline from MRA acquisition to MRA simulations[END_REF]: a validation process of the entire VIVABRAIN pipeline was performed on a physical phantom (a double bifurcation fluid circuit), on which morphological and flow data were acquired within the project. The most up-to-date full methodological framework is described in [13]. A software pipeline is associated to this framework [Ang] and open-data is a complementary purpose of the software framework, itself developed in an open-source way. It will hopefully constitutes a versatile tool for progressing in the understanding of vascular networks, especially in the brain, and the associated imaging technologies. Furthermore, it provides the foundation of a reliable and efficient computational framework allowing for a computer-aided interpretation of clinical data, as available in the Eye2brain project that will be discussed in the next sections.

Conclusions and outlook

Our contribution with respect to existing works on three-dimensional models for blood flow in realistic geometries considered the importance of boundary conditions, with a focus on the cerebral venous network. We developed the first three-dimensional sound mathematical and computational model for the venous part of the cerebral vasculature, that complements previous works on the arterial side and on reduced models for the same compartment. The framework incorporates a new Lagrange multiplier-based numerical method to take into account boundary conditions involving pressure, that we throughly analyzed in the case of the Stokes problem. We also made first steps to take into account uncertainties in different modeling assumptions through a sensitivity analysis approach, in a high performance computing framework that was validated on significant benchmarks. Finally, all these contributions were integrated in a general software pipeline for generating efficiently virtual MRA images naturally associated to a real MRA.

CHAPTER 1. THREE DIMENSIONAL MODELS

In addition to the specific extensions discussed in each section, we present hereafter some overall perspectives on the topic. Eventually, this mathematical and computational modeling framework could contribute to elucidate the numerous open questions connecting brain venous hemodynamics and neurological diseases, as recently reviewed in [START_REF] Toro | Brain venous haemodynamics, neurological diseases and mathematical modelling. A review[END_REF].

Multiscale approach. By means of the new algorithm we developed in [START_REF] Carichino | Energy-based operator splitting approach for the time discretization of coupled systems of partial and ordinary differential equations for fluid flows: the Stokes case[END_REF] (see Section 2.3), we aim at building a multiscale model that couples the 3D intracranial description of the venous network [11] with the reduced 0d model of the extra-cranial drainage pathways proposed in [START_REF] Gadda | A new hemodynamic model for the study of cerebral venous outflow[END_REF]. We made a first step by coupling the three-dimensional domain with a very simple Windkessel model in Section 1.3, that already gave promising results. However, in our current model, only two main drainage veins in the supine position are included (corresponding to the jugular veins), due to the limitations of the image acquisition process. The advantage of this new multiscale model is that it would allow to add a more sophisticated description of the collateral pathways 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] Gisolf | Human cerebral venous outflow pathway depends on posture and central venous pressure[END_REF][START_REF] Stoquart-Elsankari | A phase-contrast MRI study of physiologic cerebral venous flow[END_REF]. From the mathematical viewpoint, that would require extending our algorithm to the Navier-Stokes system, as explained in Section 2.3, and improving linear solvers efficacy by means of well-suited block-preconditioning strategies in the spirit of [START_REF] Elman | Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics[END_REF]. More generally, incorporating the possibility to describe by a reduced model a part of the entire vascular tree would be meaningful to tackle the issue of parameter identification, in line with recent works [START_REF] Lombardi | Inverse problems in 1D hemodynamics on systemic networks: A sequential approach[END_REF][START_REF] Caiazzo | Assessment of reduced-order unscented Kalman filter for parameter identification in 1-dimensional blood flow models using experimental data[END_REF], that use a 1d model of the arterial network and a the reduced-order unscented Kalman filter on synthetic or experimental data. For a further discussion on these issues see Section 2.4.

Multiphysics modeling. The proposed model for the cerebral venous system could be enhanced by coupling the current three-dimensional fluid dynamics description of the flow with new models at the inflow and the outflow. On one hand, at the inflow, blood comes from the microcirculation, comprising a huge number of small vessels. An homogenization of this vascular bed based on the Darcy law can be used in this case, but the coupling between either Navier-Stokes or a reduced model and the Darcy problem requires non-standard interface conditions and treatment, see for an extended discussion [START_REF] Angelo | On the coupling of 1D and 3D diffusion-reaction equations: applications to tissue perfusion problems[END_REF]. This approach would pave the way for the construction of a model for the whole cerebral circulation, allowing the coupling of the venous tree with the arterial network and the capillary bed. On the other hand, as discussed in Section 1.1, a rigid wall model was adopted everywhere in the network in a first approximation. However, at the outflow, the jugular vessels are more compliant that the intracranial ones, therefore a coupled fluid-structure model would be more appropriate to describe vessel-blood dynamics in this part of the network. This would require the derivation of a new model coupling the 3d Navier-Stokes system (1.1)-(1.2) on a first portion of the domain with another 3d compartment, where a fluid-structure system (see Section 3.1.1) describes the flow. To the best of our knowledge, this problem is open: we would first need to derive appropriate interface conditions to ensure the well-posedness of the problem and subsequently find an appropriate numerical strategy to solve it without the prohibitive cost of an overall fluid-structure resolution.

Mathematical and computational models for the cerebrospinal fluid (CSF). Several recent reviews [Kur11; Lin+16; BHM16] pointed out the importance of the coupling between the cerebrovascular blood flow dynamics and cerebrospinal fluid mechanics. In the framework of the VIVABRAIN project [Viv], S. Garnotel's thesis [START_REF] Garnotel | Modélisation numérique de la pression intracrânienne via les écoulements du sang et du liquide cérébrospinal mesurés par IRM de flux[END_REF] also explored different methods for the numerical modeling of the intracranial pressure and comparison with real measurements by flow MRI and intracranial pressure monitoring. Nevertheless, the understanding of the interplay be-tween blood pressure, intracranial pressure and CSF pressure (as measured by a lumbar puncture) is still scarce. To have a better description of available reduced and distributed models for CSF modeling, we are currently writing an invited book chapter entitled Mathematical modeling of the cerebrospinal fluid flow, in collaboration with L. Sala (PhD student, Univ. de Strasbourg) and F. Salerni (PhD student, Univ. of Parma), that reviews different contributions on this topic. The chapter is included in the book Mathematical modeling of ocular fluid dynamics: From theory to clinical applications (A. Harris and G. Guidoboni Eds., Springer) and gives a literature perspective to the recent eye-cerebral model we developed in [8; 17], see also Section 2.1.

Chapter 2

Reduced and multiscale mathematical and computational models for biofluids. Applications to the study of the coupled eye-cerebral system Biofluid models and simulations in patient specific geometries and under physiological conditions at the macro-scale are usually based on a three-dimensional description of the flow by means of a system of partial differential equations (PDEs). To simplify the numerical simulations and the interpretation of the results, various modeling reduction techniques have been proposed, see for instance [FQV09, Chap.10] and the references therein. They are able to capture the main physical phenomena governing the system, at a lower computational cost than the full solution of a PDEs problem.

One-dimensional models take advantage of the basically cylindrical morphology of biological conduits. They can be derived by resorting to different procedures, such as cross-sectional averaging techniques or asymptotic analysis methods. The resulting mathematical description is based on 1d hyperbolic PDEs, that can capture wave-propagation phenomena in the entire circulatory system with a reasonable computational cost. This is not the approach we will follow in the sequel, since we rather focused on an even more simplified, 0d description of the flow.

Reduced 0d models provide a circuit-based representation of the fluid dynamics in each compartment, based on the analogy between electric and hydraulic networks, see Table 2.1. By writing Kirchhoff laws for the nodes (conservation of current/flow rate) and for closed circuits (conservation of the voltage/pressure), the resulting mathematical model is a system of differential algebraic equations, potentially nonlinear. These models, also called lumped parameter models, are used in different ways depending on the specific modeling needs, as recently reviewed in [START_REF] Quarteroni | Geometric multiscale modeling of the cardiovascular system, between theory and practice[END_REF]. In a first approach, many studies implemented them as boundary conditions for three-dimensional blood flow simulations in regions of particular interest. Alternatively, 0d reduced models have been used to provide systemic descriptions of closed circulatory systems, where three-dimensional regions might be embedded. From this so-called geometric multiscale modeling perspective, the PDE-ODE coupling leads to challenging theoretical and numerical issues, since appropriate interface conditions enforcing the continuity of mass and the balance of forces need to be devised and should also be preserved at the discrete level. In this research direction, the purpose of the chapter is two-fold: (i) first, we present in Sections 2.1 and 2.2 our contributions towards the elaboration of several reduced 0d models describing the coupled dynamics of different biofluids in the eye-cerebral system; (ii) second, we describe in Section 2.3 a new methodological strategy for the numerical solving of coupled systems of partial and ordinary differential equations for fluid flows. The latter was motivated by applications to blood flow modeling through the cardiovascular system from a multiscale perspective, even though the resulting conceptual framework may be meaningful and applicable to a more general context of hydraulic networks.

Design and development of a lumped model for the coupled eye-cerebral system

In the previous chapter, progress made in understanding cerebral hemodynamics by means of mathematical and computational models highlighted an increasing need for clinical data to be used in the calibration and validation process. To mitigate measurement problems in the brain, we are considering here an alternative approach: the eye is the only place in the human body where structural and functional vascular features can be observed and measured easily and noninvasively down to the capillary level. Many clinical studies have shown correlations between hemodynamic alterations in the eye and pathological conditions in the brain [START_REF] Heringa | Associations between retinal microvascular changes and dementia, cognitive functioning, and brain imaging abnormalities: a systematic review[END_REF][START_REF]Interconnection Between Brain and Retinal Neurodegenerations[END_REF]. This is not surprising, since the two organs share anatomic, embryologic, and physiologic characteristics [START_REF] Heringa | Associations between retinal microvascular changes and dementia, cognitive functioning, and brain imaging abnormalities: a systematic review[END_REF], see Figure 2.1. However, the mechanisms underlying these correlations are still elusive, likely due to the numerous factors that influence the relationship between eye and brain and the intrinsic difficulty of isolating these factors in a clinical setting and measuring their individual contribution. These factors include arterial blood pressure (BP), intraocular pressure (IOP), intracranial pressure (ICP) and cerebrospinal fluid pressure (CSF-p).

The goal of this section is to propose a new mathematical model describing the coupled dynamics of different biofluids in the brain and in the eye. To the best of our knowledge, an approach describing at this level of details the interactions between all the compartments described below was not available in the literature. Details on the design and development of the model are presented in Section 2.1.1. In Section 2.1.2, the model is employed to investigate how microgravity conditions affect ocular and cerebral fluid-dynamics. The presentation is based on [8; 17] and a full paper in preparation. It also benefits from the progress we made in understanding the aqueous humor dynamics through mathematical and computational models, as reviewed in Section 2.2.

Construction of the model

We illustrated in Figure 2.1 the complex interplay between several biofluids that are involved in the fluid-dynamic connections between the eye and the brain. In a first modeling attempt, we proposed a reduced 0d description of the eye-cerebral system, connected with a very simplified model of the body, the electrical analogue of which is shown in Figure 2.2. The physiological network is subdivided into a number of linked, interacting compartments, each of which contains a single physical constituent, such as blood, cerebrospinal fluid, interstitial fluid and aqueous humor.

The interactions among the following components are incorporated (see Figure 2.2):

in the brain:

• blood (3 compartments: intracranial arteries (I), capillaries (C), intracranial venous sinuses (S));

• cerebrospinal fluid (1 compartment: ventricular CSF (F));

• cerebral tissue and interstitial fluid (1 compartment: brain (B));

in each of the two eyes (see Figure 2.3):

• blood (5 compartments: retina (r ), choroid (ch), central retinal artery (cra), central retinal vein (crv ), ciliary circulation (cl ));

• aqueous humor (ah) (1 compartment: anterior and posterior chamber); and finally

• cerebrospinal fluid (1 compartment: extra-ventricular (T) bridging intracranial and extracranial regions, also including the subarachnoid space in the optic nerve posterior to the lamina cribrosa).

P A R A P V R V Q CF R I R I R C R C R S R S R F T R T S R BT R F B R T V K CB Ventricular CSF (F) Extra-ventricular CSF (T) Brain Tissue (B) Intracranial Arteries (I) Capillaries (C) Venous Sinuses (S) Q CB Q F B Q BT Q CS Q IC Q T S Q F T Q OA Q Eye,In Q Eye,In Q CvS Q Eye,Out Q Eye,Out Ophthalmic Artery R I,OA Cavernous Sinus BRAIN Arteries (A) Veins (V) BODY Lamina cribrosa ciliary (cl) choroid (ch) retina (r) (cra) (crv) Aqueous humour (ah) EYE Lamina cribrosa ciliary (cl) choroid (ch) retina (r) (cra) (crv) Aqueous humour (ah) EYE Q A Q T V Q SV Q V Figure 2.
,v R ch,in R ch,a1 R ch,a2 R ch,c1 R ch,c2 R ch,v1 R ch,v2 R ch,out Pch,out Pcl,in Pcl,a Pcl,c Pcl,v R cl,in R cl,a1 R cl,a2 R cl,c1 R cl,c2 R cl,v1 R cl,v2 R cl,out Pcl,out Lin Jsecr IOP Rtm Episcleral Veins Pevp Ruv 1 2 3 4 5 QEye,in Q ch Q cl Qr QEye,out Q ch Q cl,v Qr J uf J ah Juv Jtm

CHAPTER 2. REDUCED AND MULTISCALE MODELS

To reflect the fluid dynamic connections between the eye and the brain, the models are coupled through the following relations:

• the node corresponding to the ophthalmic artery (OA) pressure in the eye model (denoted P Eye,in in Figure 2.2) is connected to the intracranial arterial pressure in the brain model (P I ) through an effective vascular resistance (R I,OA );

• the node corresponding to the cavernous sinus pressure (P Eye,out ) and that corresponding to the episcleral venous pressure (P evp ) in the eye are directly connected to the venous sinus pressure in the brain model (P S );

• the lamina cribrosa is acted upon by the IOP from the ocular side (light blue arrow) and by the CSF pressure in the subarachnoid space in the optic nerve, here assumed to be equal to the extraventricular CSF pressure (grey arrow).

The flow is driven by the following mechanisms:

• the pressure drop between the central arteries (P A ) and the central veins (P V );

• the active secretion of aqueous humor in the eye (J secr );

• the production of CSF in the brain (Q CF ).

Remark 2.1.1. In this work we focus on steady simulations, neglecting, in particular, time variations occurring on the time scale of heart beat. We also neglect autoregulation mechanisms of small vessels, since we wish to keep the model relatively simple in order to understand its basic behavior.

By writing the Kirchoff law of currents at the circuit nodes, we obtain a set of nonlinear algebraic equations. The nonlinearity is a consequence of the fact that, in some compartments, resistances are assumed to depend on the pressures. Following [START_REF] Guidoboni | Intraocular Pressure, Blood Pressure, and Retinal Blood Flow Autoregulation: A Mathematical Model to Clarify Their Relationship and Clinical Relevance[END_REF] and [START_REF] Pedley | The Fluid Mechanics of Large Blood Vessels[END_REF], we model these deformable tubes as Starling resistors, reflecting the physiological high collapsibility of these vessels when the transmural pressure becomes negative. All constitutive laws are listed in Table 2.2, and two main types of flows are included in the model: filtration and pressure-driven flows. See for more details on the mathematical modeling approach Section 2.2.1. Precise values for model parameters, the full set of notations and the validation procedure are available in the complete version of the paper.

We are not aware of works in which a mathematical model describing the complex relationship between different biofluids in the brain and in the eye in such a detailed manner has been proposed. The lumped parameter circuit for the brain stems from [START_REF] Lakin | Modelling the Response of Intracranial Pressure to Microgravity Environments[END_REF], where a mathematical model was employed for microgravity simulations. The eye model originates from [START_REF] Guidoboni | Intraocular Pressure, Blood Pressure, and Retinal Blood Flow Autoregulation: A Mathematical Model to Clarify Their Relationship and Clinical Relevance[END_REF] for the retinal circulation and has been extended to account for the three ocular vascular beds (retina, choroid and ciliary body) adapting work from [START_REF] Kiel | Ciliary blood flow and aqueous humor production[END_REF]. Finally, ocular hemodynamics is coupled with the aqueous humor dynamics, as modeled in [START_REF] Szopos | Mathematical modeling of aqueous humor flow and intraocular pressure under uncertainty: towards individualized glaucoma management[END_REF] and extended in [14]. The advantage of this new network-based model is that it provides a systemic view, able to capture the overall dynamics of the interwoven physiology of blood, CSF, ocular humors and interstitial fluids in the eye and in the brain, while maintaining a relatively accessible mathematical complexity and low computational costs.

Microgravity conditions modeling and numerical results

Understanding how microgravity conditions affect ocular and cerebral fluid-dynamics and the role that the vascular components may play in the (potentially permanent) loss of visual acuity after space flight is a challenging issue. This problem is called Visually Impairment and Intracranial Pressure (VIIP) syndrome [START_REF] Nelson | Microgravity-Induced Fluid Shift and Ophthalmic Changes[END_REF]. As part of several clinical and research projects directed in the International Space Station, NASA recently initiated an extended acquisition campain to analyze and identify astronaut-specific anatomical and physical factors that may contribute to increase the risk of developing VIIP during long-duration space flights. Many factors are hypothesized to contribute to VIIP, among which upper body fluid shift, alterations in intraocular pressure, intracranial pressure, tissue biomechanics and blood flow. Due to the difficulty to single out each of these factors using in-vivo studies, we proposed a complementary approach, based on the novel mathematical model described in the previous section, to evaluate the interactions between fluid flows and pressures in the brain and eyes and their mechanical implications in VIIP.

To our knowledge little work has been done concerning the mathematical modeling of VIIP.

The work [START_REF] Nelson | Computational Modeling of Cephalad Fluid Shift for Application to Microgravity-Induced Visual Impairment[END_REF] presented a lumped-parameter model to calculate mean ICP, IOP and blood flows in the eye in zero gravity conditions. However, in the authors' model, the eye is treated in a simplified manner, as a distensible compartment inside of which a second distensible compartment exists, representing the intraocular vasculature. The contribution [START_REF] Vera | Lumped Parameter Models of the Central Nervous System for VIIP Research[END_REF] developed a lumpedparameter model of fluid transport in the central nervous system aimed at simulating the influence of microgravity on ICP. As a forthcoming step, the authors' intention is to couple this model to lumped parameter and finite element models of the eye.

We performed simulations using the model described in Section 2.1.1 in four conditions, considered to be relevant in VIIP [START_REF] Lakin | Modelling the Response of Intracranial Pressure to Microgravity Environments[END_REF], for different imposed blood/aqueous and blood/brain osmotic pressure differences Δπ:

• LHDT (long term head down tilt, a ground-based experimental procedure that is used to simulate the effects of microgravity on the cardiovascular system): we account for a small gravitational effect, which increases the central venous and arterial pressures;

• M0 (microgravity with an intact blood/brain barrier): we set the gravitational contribution and the central venous pressure equal to zero;

• M1, M2 (microgravity with a weakened blood/brain barrier): in addition to the previous case, we also increase the permeability of the blood/brain barrier by a factor 1.75 in the case M1 and by a factor 2 in the case M2.

As an illustration, we present in Figure 2.4 the results of simulations for IOP, ICP and compressive stress in the lamina in the four simulated cases; they show that a decrease in Δπ leads to: (i) an approximately linear increase in ICP, more marked if the blood/brain barrier is weakened, as in [START_REF] Lakin | Modelling the Response of Intracranial Pressure to Microgravity Environments[END_REF]; (ii) a nonlinear increase in IOP, due to the venous collapse as the transmural pressure becomes negative. Simulated ocular blood flow Figure 2.5 shows that, as Δπ decreases: (i) ciliary and choroidal blood flows decrease almost linearly, until vessel collapse occurs, and then do so much more rapidly; (ii) choroidal circulation the flux in the retina initially grows owing to a purely perfusion control mainly enacted by the intraocular venous segments.

Conclusions and perspectives.

From the clinical viewpoint, these results showed interesting features. Specifically, when simulating microgravity conditions, the model predicted that: (i) the

CHAPTER 2. REDUCED AND MULTISCALE MODELS

venous segments play a fundamental role in controlling pressures and fluxes in the ocular circulation, owing to the possibility of their collapse; (ii) IOP and ICP increase in microgravity conditions, but their respective trends may be different; (iii) retinal circulation is less susceptible to microgravityinduced alterations than that in choroid and ciliary. These findings aid VIIP understanding, indicating that clinical assessment of ocular venous function may be considered as a determinant factor, for which astronauts could be screened on earth and in-flight. We finally remark that the ingredients embedded in the proposed model make it a very flexible tool to study interactions between the eyes and the rest of the body (particularly the brain). Thus, besides the application to VIIP, we believe that the model could serve as a framework on which more can be built and could be easily extended to study other pathological states, such as glaucoma (see also Section 2.2.1).

From the mathematical modeling standpoint, the outcome of this study is the first mathematical description of the biophysical connection between the eye and the brain. The model was implemented in the open-source modeling and simulation environment OpenModelica [Ope] and validated against clinical data. Next steps include (i) a sensitivity analysis study to theoretically identify systemic, ocular and cerebral parameters, in the spirit of our recent contributions [14; 21] and (ii) the extension of the model to the time-dependent case. The major challenge of the latter is the strongly nonlinear multiscale temporal nature of the problem, that in addition to the nonlinearities coming from the constitutive laws will call for well suited numerical methods for highly stiff nonlinear ODE systems.

A zoom on mathematical models for the aqueous humor

We present in this section our work on mathematical models for the aqueous humor dynamics and their coupling with ocular hemodynamics, as part of the complex eye-cerebral model described in Sec. 2.1. The aqueous humor is a transparent, gelatinous fluid similar to plasma located in the anterior and posterior chambers of the eye. The balance between production and drainage of aqueous humor regulates the fluid pressure inside the eye, also known as intraocular pressure (IOP). Elevated IOP is clinically referred to as ocular hypertension and it represents a major risk factor for irreversible vision loss, as in glaucoma. However, the establishment of optimal IOP levels for patients is still controversial, as recently reviewed in [START_REF] Vercellin | Mathematical modeling and glaucoma: the need for an individualized approach to risk assessment[END_REF]: clinical studies report evidence of disease progression despite a significantly reduced IOP level and optical nerve damage and vision loss despite normal IOP level. Therefore, understanding what is the relative contribution of each risk factor and/or different combination of risk factors is of major importance. Another meaningful open question is related to what determines the efficacy of IOP-lowering medications. To date, this is the only modifiable risk factor [START_REF] Vercellin | Mathematical modeling and glaucoma: the need for an individualized approach to risk assessment[END_REF]. However, efficacy varies among individuals,e.g: the lowering effect is more significant when IOP is higher and it changes depending on the time of day, but the underlying mechanisms are not yet fully understood.

Motivated by these considerations, we present in Section 2.2.1 a mathematical model that describes the steady state value of IOP under uncertainty. We subsequently performed numerical simulations based on this model to assess the outcome of IOP-lowering medication used in glaucoma management. In Section 2.2.2 we further propose two increasingly complex models to couple aqueous humor production and drainage with ocular blood circulation.

Modeling the intraocular pressure under uncertainty

Only a few modeling works have studied aqueous humor flow and its relation to IOP-lowering medications, e.g. [START_REF] Lyubimov | Dynamics of the Intraocular Fluid: Mathematical Model and Its Main Consequences[END_REF][START_REF] Kiel | Ciliary blood flow and aqueous humor production[END_REF]. Importantly, none of them explicitly accounted for uncertainties and variabilities in the model parameters. We proposed in [START_REF] Szopos | Mathematical modeling of aqueous humor flow and intraocular pressure under uncertainty: towards individualized glaucoma management[END_REF] a simplified mathematical model describing the balance between aqueous humor production and drainage; we then performed a sensitivity analysis aimed at quantifying the influence of parameters' variations on the IOP distribution in various situations of clinical interest. Preliminary findings on this topic were published as peerreviewed conference abstracts in [7; 2].

Description of the model. The mathematical approach is based on the analogy between transport of mass and flow of an electric current. We give hereafter the main steps in the derivation of the model, as an illustrative simple example of our approach when constructing lumped-parameter models for this type of biological problems and how different phenomena were incorporated. The network-based view is schematized in Figure 2.

(blue portion).

Aqueous humor is produced at the level of the ciliary body by a combination of a passive mechanism, the ultrafiltration, and an active mechanism, the ionic secretion, and is modulated by the total inflow facility L [LMS07; Kie+11]. Here the term facility indicates hydraulic conductance, namely a flow rate per units of pressure. The ultrafiltration from the ciliary circulation consists of flow of transparent fluid across semipermeable membranes and is driven by blood/aqueous humor differences in hydrostatic pressures (cBP -IOP ) and oncotic pressures Δπ p : the latter is modulated by a protein reflection coefficient σ p . The notation cBP stands for the blood pressure in the capillaries of the ciliary body. The inflow, as a result of the active ionic secretion, is proportional to the blood/aqueous humor osmotic pressure difference Δπ s , via a reflection coefficient for low-molecular components σ s .

The drainage of aqueous humor from the eye is driven by passive mechanisms through two different pathways. The trabecular pathway, also known as conventional pathway, consists of aqueous humor flow through the trabecular meshwork, into the Schlemms canal and the episcleral veins. As proposed in [START_REF] Brubaker | The effect of intraocular pressure on conventional outflow resistance in the enucleated human eye[END_REF], the model consists of a flow through a nonlinear resistor positioned between the anterior chamber (where pressure is equal to IOP) and the episcleral veins (where pressure is equal to EVP), with outflow facility depending on R 0 (resistance when IOP equals EVP) and Q, an outflow obstruction coefficient. The uveoscleral pathway, also known as the non-conventional pathway, consists of aqueous humor flow through the ciliary muscle and into the supraciliary space. This contribution is modeled as the flow through a nonlinear resistor connected to the ground [START_REF] Kiel | Ciliary blood flow and aqueous humor production[END_REF], with an outflow facility depending nonlinearly on the pressure through the Michaelis-Menten type relation [START_REF] Johnson | The original Michaelis constant: translation of the 1913 Michaelis-Menten paper[END_REF]. The constant k 1 is the maximum value attainable by the uveoscleral flow rate and k 2 is the Michaelis constant for the uveoscleral flow rate, namely the pressure value for which the uveoscleral flow rate is half of k 1 .

Consequently, the steady state value of IOP, resulting from the balance between production and drainage of aqueous humor can be written as:

L[(cBP -IOP )-σ p Δπ p -σ s Δπ s ] = 1 R 0 [1 + Q(IOP -EV P )] (IOP -EV P )+ k 1 k 2 + IOP IOP. (2.1)
This is a scalar third-order polynomial equation in the sole unknown and can be explicitly computed from the previous formula. Control state values for the parameters and more details on the CHAPTER 2. REDUCED AND MULTISCALE MODELS derivation of the model can be found in [START_REF] Szopos | Mathematical modeling of aqueous humor flow and intraocular pressure under uncertainty: towards individualized glaucoma management[END_REF].

To include potential sources of uncertainties as well as to identify and rank parameters having the most important influence on IOP, we applied a global stochastic sensitivity analysis to the model described above. We considered stochastic variations in cBP following a normal distribution, and in L, Δπ s , C 0 = 1 R 0 (trabecular outflow facility), k 1 and EVP following a uniform distribution, both within physiological ranges. By using the probability distribution of IOP, we computed variance-based sensitivity indexes, also known as Sobol indexes [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF] and the probability density function [START_REF] Saporta | Probabilités, analyse des données et statistique[END_REF] which describes the relative frequency of a given IOP value. For each parameter, its direct influence on is quantified in terms of first-order Sobol indexes, and the influence through interactions with other parameters is identified by means of the total Sobol indexes. The values of first-order and total indexes can be estimated via Monte Carlo simulations [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF] or via reduced order models using polynomial chaos expansion [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF]. The former method is very costly from the computational viewpoint as it requires many evaluations to ensure convergence, whereas the latter requires considerably less evaluations. Both methods have been compared and provide similar results. We report in the sequel the results obtained using the polynomial chaos reduced model. Computations were performed using Matlab/Octave and the open-source software platform OpenTURNS [START_REF] Baudin | OpenTURNS: An Industrial Software for Uncertainty Quantification in Simulation[END_REF].

Numerical results and discussion. The model was used to compute the distribution in four different cases of clinical interest: (i) ocular normotensive healthy subjects (ONT); (ii) ocular hypertensive subjects (OHT); (iii) ONT subjects treated with IOP-lowering medications (ONTm); and (iv) OHT subjects treated with IOP-lowering medications (OHTm).

We only report here the results for the first category (healthy individuals), that also validate the model by comparison with experimental results, see Figure 2.6. The probability density function (left panel) fits a right-skewed Gaussian curve with a frequency peak of 25% at 15.13 mmHg and a skewness of 0.2, which is in a very good agreement with the results from a population-based study on approximately subjects 12.000 subjects [START_REF] Carel | Association between ocular pressure and certain health parameters[END_REF] (green curve in Figure 2.6, left panel). The results for the Sobol indexes (Figure 2.6, right panel) suggest that IOP is strongly influenced by cBP and and Δπ s and mildly influenced by the levels of L, C 0 and EV P . The influence of k 1 on IOP appears to be minimal.

The outcomes of the simulations in the other three cases are discussed in detail in [START_REF] Szopos | Mathematical modeling of aqueous humor flow and intraocular pressure under uncertainty: towards individualized glaucoma management[END_REF]. To summarize, the proposed model suggests that the outcomes of IOP-lowering treatments depend on the initial IOP level of the patient and on its individual clinical condition. Our analysis also suggests that IOP-lowering effects are more pronounced when aqueous humor production is affected rather than aqueous humor drainage. The effects of lowering IOP are also more apparent when trabecular outflow is increased instead of the uveoscleral outflow. Another interesting finding of our analysis is that a patients blood pressure strongly influences the outcomes of IOP-lowering treatments, which may explain why the effect of some drugs differ between day-time and night-time and/or amongst individuals.

Conclusions and perspectives. Elevated IOP is a recognized risk factor for vision loss. Numerous risk factors are suggested in the literature [START_REF] Vercellin | Mathematical modeling and glaucoma: the need for an individualized approach to risk assessment[END_REF]: age, gender, ethnicity, genetic factors, vascular dysregulation, diabetes, blood pressure etc. However, the establishment of optimal IOP levels for patients is still controversial and what is the relative contribution of each risk factor and/or different combination of risk factors is still not well understood. We proposed a mathematical model to theoretically investigate the relationship between several uncertain parameters influencing the flow and evaluated their potential impact on the individual response to IOP-lowering medications.

A further investigation that incorporates a theoretical and numerical model coupling aqueous humor production and drainage with ocular blood flow is presented in the next section. Changes in ocular blood volume, mainly localized in the choroid, are also conjectured to affect the time variations of IOP, but they were not considered in [START_REF] Szopos | Mathematical modeling of aqueous humor flow and intraocular pressure under uncertainty: towards individualized glaucoma management[END_REF]. In this direction, we recently made a first attempt to consistently combine the static and dynamic components contributing to the pressure distribution inside the eye in [START_REF] Stefanoni | Clinical assessment of intraocular pressure: a whole eye dynamic model[END_REF], by incorporating in the model a nonstationary component, mainly due to blood flow oscillations. A more in-depth analysis of the mathematical properties of the steady and unsteady versions of the model is another interesting question, since it could bring new insights on the qualitative behavior of the solutions.

Modeling the coupled dynamics of ocular blood flow and aqueous humor

As anticipated in the previous section, we present hereafter two increasingly complex mathematical models (M1 and M2) describing the coupled dynamics of aqueous humor and ocular blood flow.

We performed a sensitivity analysis in both cases with the goal of investigating the influence of different factors on the intraocular pressure production and drainage. The presentation is based on [14].

Description of the model. As previously, we adopted here a reduced-modeling approach, stemming from the analogy between electric and hydraulic networks. Model M1 couples a simplified version of the model [START_REF] Szopos | Mathematical modeling of aqueous humor flow and intraocular pressure under uncertainty: towards individualized glaucoma management[END_REF] (see Section 2.2.1 and [LMS07]) for aqueous humor production and drainage with the models for the blood flow in the ciliary, choroidal and retinal circulations proposed in [START_REF] Kiel | Ciliary blood flow and aqueous humor production[END_REF] and [START_REF] Guidoboni | Intraocular Pressure, Blood Pressure, and Retinal Blood Flow Autoregulation: A Mathematical Model to Clarify Their Relationship and Clinical Relevance[END_REF].

In model M2, the blood flow description is the same as in M1, whereas the aqueous humor dynamics model includes more complex phenomena, following [START_REF] Kiel | Ciliary blood flow and aqueous humor production[END_REF]: plasma filtration and biochemical reactions are considered in the multi-layered region of the eye comprising the stroma and the epithelial ciliary layer, leading to the posterior chamber. The model is schematized in Fig. 2.7.

CHAPTER 2. REDUCED AND MULTISCALE MODELS

Enforcing the Kirchhoff Current Law at each node of the circuit, except where given data are assumed and where connection to ground is applied, one obtains a nonlinear system of algebraic equations in the unknown vector where C st and C aq are the protein concentrations in the stroma and in the aqueous humor, P cil c is the blood pressure in the ciliary body capillaries, P st , P post and P ant are the aqueous humor pressures in the stroma and in the posterior and anterior chambers, m bl basolateral ionic concentration, and P O2 c , P O2 v , P O2 ce , P O2 aq are the partial pressures of oxygen in the ciliary body capillaries, veins, epithelium and aqueous humor. Thus, the overall models M1 or M2 consists in two main blocks, one for the aqueous humor flow and one for blood flow, which are strongly coupled, since: (i) the blood pressure in capillaries and veins contributes to regulate production and drainage of aqueous humor, and (ii) the aqueous humor pressure acts as an external pressure on the blood vessels thereby altering blood flow. For the solution of the overall problem a block fixed point iteration is adopted until convergence.

P = � C st ,
� � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � ���� ���� ���� ���� ����� ����� ����� ����� ���� ���� ���� ���� ���� ���� � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � G 1 G 2 I 1 I prot out G c G ce
Numerical results. Computations of the first order Sobol indexes [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF] based on M1 confirmed, as in [START_REF] Szopos | Mathematical modeling of aqueous humor flow and intraocular pressure under uncertainty: towards individualized glaucoma management[END_REF], that acting on the aqueous humor inflow seems more effective than changing aqueous humor outflow. In addition, regarding hemodynamic aspects, simulations showed that pressure in ciliary capillaries and retinal blood flow are strongly influenced by blood pressure and systemic vascular resistance, which are systemic factors depending on the general condition of the individual.

We also performed a sensitivity analysis on model M2, but only in a deterministic framework. Because of the increased complexity of the formulation, a self-consistent coupling with the blood flow model was not addressed in a preliminary validation, and computations were performed by assuming the dependent variables P cil c , P O2 c and P O2 v in (2.2) to be given quantities. In particular, we focused our attention in the computation of P st and P ant (P post is assumed to be equal to P ant ), and of the pressure P vitreous in the vitreous which coincides with the intraocular pressure. Simulations were run to evaluate the dependence of these three variables on the variation (i) of the production of aqueous humor P AH and (ii) of the total aqueous humor drainage D AH . A complete description of the parameters and results is available in [14]. As an illustration, we report in Figure 2.8 the outcomes of the simulations in the case where the drainage of aqueous humor D AH is varied. Results show that the intraocular pressure is less sensitive to parameter variation, unlike stromal and anterior chamber pressures. In particular, we see that, as physically expected, these latter pressures decrease appreciably as D AH becomes larger. However, results clearly reveal that a further increase of drainage has no positive effect on pressure decay. This is reflected into a saturation of all the pressure curves to a limit value, as visible in the right panel of Figure 2.8.

Discussion and perspectives. The mathematical models M1 and M2, coupling aqueous humor production and drainage with ocular blood flows, lead to a better understanding of the delicate, yet important, relationship between different factors determining intraocular pressure. Further on, this theoretical investigation helped in identifying patients-specific factors (as for instance systemic blood pressure) that influence the efficacy of different medications used in glaucoma management. The complex multi-physical nature of this coupled dynamics called in a first approximation for several simplifications: the model is steady and does not account for autoregulation mechanisms in small vessels. In this direction, the aqueous humor block could incorporate a nonstationary component as we recently proposed in [START_REF] Stefanoni | Clinical assessment of intraocular pressure: a whole eye dynamic model[END_REF]; as for the ocular blood flow contribution, time-dependent models including autoregulation in the retina, with different levels of complexity, are developed in [Gui+14; AGL16] and could be coupled with the first block. From the mathematical viewpoint, taking into account the different nonlinearities due to autoregulation phenomena raises new interesting questions both at the theoretical and computational level; another challenge comes from the understanding of the unsteady component of the system and its asymptotic behavior.

A different development was triggered by recent discussions with A. Mauri and R. Sacco, that aims at coupling the current model that describes aqueous humor dynamics at the ocular tissue scale model (coarse level) to an epithelial cellular scale model (fine level). More precisely, the idea is to keep the description of the drainage path and of the ultrafiltration process at a macroscopic level and to expand the model of the ionic pump at the fine level in the spirit of the recent model proposed by my colleagues and their collaborators in [START_REF] Mauri | Electro-fluid dynamics of aqueous humor production: simulations and new directions[END_REF]. Such an approach would hopefully provide insights in understanding the role of the sodium-potassium pump as an important player in the aqueous humor production, not only at the cellular level, but also at a macroscopic scale. In [START_REF] Mauri | Electro-fluid dynamics of aqueous humor production: simulations and new directions[END_REF], the velocity-extended Poisson-Nernst-Planck equations describing ion electrodiffusion are coupled to the Stokes equations to describe aqueous humor flow into the basolateral space adjacent to the nonpigmented ephitelial cells. The coupling between this three-dimensional partial differential equations system and the reduced 0D description of the rest of the circuit presents interesting modeling (for boundary and interface conditions) and numerical (for the solving algorithm) issues. In particular, the new numerical approach based on a splitting procedure that we detail in the next section could provide a sound starting point in devising an appropriate algorithm for the problem.

A new algorithm for the numerical solving of coupled distributed-0d problems

Many strategies have been adopted for the numerical solution of coupled PDE-ODE systems in the context of blood flow modeling. In particular, monolithic and splitting (or partitioned) schemes have been proposed, where the PDE and ODE systems are solved simultaneously or in separate substeps, respectively. An extensive discussion about their advantages and limitations can be found in [START_REF] Quarteroni | Geometric multiscale modeling of the cardiovascular system, between theory and practice[END_REF].

The contribution presented in this section focuses on splitting techniques and on the properties of their modular structure. Among the many interesting contributions in this area, we mention here those that are most closely related to our work. In [START_REF] Quarteroni | Coupling between lumped and distributed models for blood flow problems[END_REF], the multiscale coupling between the Navier-Stokes equations in a rigid domain and a lumped parameter model is considered. A splitting strategy based on subiterations between PDE and ODE solvers at each time step is proposed and assessed in different meaningful configurations. This splitting formulation was then used as an effective tool to prove the well-posedness of the coupled problem [START_REF] Quarteroni | Analysis of a geometrical multiscale model based on the coupling of ODE and PDE for blood flow simulations[END_REF], in combination with appropriate fixed point results. In [START_REF] Fouchet-Incaux | Numerical stability of coupling schemes in the 3d/0d modelling of airflows and blood flows[END_REF], the authors compare the numerical stability of explicit and implicit coupling between the Stokes or Navier-Stokes equations and circuit-based models containing resistances and capacitances. Unconditional stability was proved in the implicit case, whereas conditional stability was proved in the explicit case. In [START_REF] Moghadam | A modular numerical method for implicit 0D/3D coupling in cardiovascular finite element simulations[END_REF], a time implicit approach to couple general lumped parameter models with a finite-element based solution of a Navier-Stokes problem in a 3d domain is proposed. The algorithm combines the stability properties of monolithic approaches with the modularity of splitting algorithms. The method is based on a Newton type iterative scheme, where data are exchanged between the two domains at each Newton iteration of the nonlinear Navier-Stokes solver to ensure convergence of both domains simultaneously.

To the best of our knowledge, the splitting schemes that have been proposed for coupled PDE-ODE systems in the context of fluid flow modeling so far, require subiterations between substeps, usually involving the values of pressure and flow rate at the multiscale interfaces, in order to achieve convergence of the overall algorithm. Depending on the mathematical properties of the models in each substep and of the coupling between them, the convergence of such subiterations might become an issue, especially in the case of nonlinear problems.

We describe in this section a new method, presented in [START_REF] Carichino | Energy-based operator splitting approach for the time discretization of coupled systems of partial and ordinary differential equations for fluid flows: the Stokes case[END_REF] and detailed in a full subsequent paper, that aims at providing a novel splitting scheme for coupled PDE-ODE systems for fluid flow that does not require subiterations between substeps in order to achieve stability for the overall algorithm. The stability properties of the scheme follow from ensuring that the physical energy balance is maintained at the discrete level via a suitable application of operator splitting techniques to semi-discretize the problem in time [Glo03, Chap. II]. As a result, the proposed algorithm allows us to: (i) solve in separate substeps potential nonlinearities within the systems of PDEs and/or ODEs; (ii) maintain some flexibility in choosing the numerical method for the solution of each subproblem; (iii) ensure unconditional stability without the need of subiterations between substeps.

Description of the coupled Stokes-0d problem

Geometrical architecture of the coupled system. We start by describing the main blocks that constitute the model under consideration. The geometrical setting of the coupled problem consists of: (i) L regions of space denoted by Ω l ⊂ R d , with l ∈ L = {1, . . . , L} and d = 2 or 3, where the fluid flow is described by the Stokes equations; (ii) M lumped hydraulic circuits denoted by Υ m , with m ∈ M = {1, . . . , M }, where the fluid flow is described by the hydraulic analog of the Kirchoff laws of currents and voltages. We assume that the boundary of each domain Ω l , denoted by ∂Ω l , is the union of three portions, namely

∂Ω l = Γ l ∪ Σ l ∪ S l (2.3)
where different types of boundary and interface conditions are imposed. Specifically, Dirichlet conditions are imposed on Γ l , Neumann conditions are imposed on Σ l and Stokes-circuit coupling conditions are imposed on S l , as described below. We assume that Dirichlet conditions are imposed on a portion of the boundary of each domain Ω l , i.e. Γ l � = ∅, and that each domain Ω l is connected to (at least) one circuit, i.e. S l � = ∅. In particular, each region Ω l may have j Ω l ≥ 1 Stokes-circuit connections, implying that each boundary portion S l may be written as

S l = � m∈M l S lm (2.4)
with l ∈ L. We remark that, for each l ∈ L, the set M l ⊆ M identifies the circuits Υ m that are connected to Ω l . Similarly, each circuit Υ m may have j Υm ≥ 1 Stokes-circuit connections and the set L m ⊆ L identifies the Stokes regions Ω l that are connected to Υ m . We will also consider cases where Neumann conditions are not imposed on the boundary of Ω l , i.e. Σ l = ∅, as it happens in Example 3 described below. It may also happen that the same Stokes region Ω l and the same CHAPTER 2. REDUCED AND MULTISCALE MODELS circuit Υ m enjoy multiple connections. Thus, an additional subscript is introduced to distinguish between the various connections, so that we can write:

S lm = j Ω l ,Υm � k=1 S lm,k , (2.5)
where j Ω l ,Υm is the total number of connections between Ω l and Υ m . Note that with these notations, we also have

j Ω l = � m∈M l j Ω l ,Υm and j Υm = � l∈Lm j Ω l ,Υm . (2.6)
Stokes problems in Ω l . Let v l = v l (x, t) and p l = p l (x, t), for l ∈ L, denote the velocity vector field and the pressure field, respectively, pertaining to the fluid in each domain Ω l × (0, T ), with Ω l ⊂ R d , d = 2, 3 and T > 0. Then, for l ∈ L, we can write the Stokes equations as:

∇ • v l = 0 in Ω l × (0, T ) (2.7) ρ ∂ v l ∂t = -∇p l + µΔv l + ρf l in Ω l × (0, T ) (2.8)
where ρ and µ are positive given constants representing the fluid density and dynamic viscosity, respectively, and f l are given body forces per unit of mass. The system is equipped with the initial conditions:

v l (x, t) = v l,0 (x) in Ω l (2.9)
and the boundary and interface conditions:

v l = 0 on Γ l × (0, T ) (2.10) � -p l I + µ∇v l � n l = -p l n l on Σ l × (0, T ) (2.11) � -p l I + µ∇v l � n lm,k = g lm,k on S lm,k × (0, T ) (2.12)
where I is the d × d identity tensor, n l is the outward unit normal vector to Σ l , and p l = p l (t) are given functions of time. For m ∈ M l and k = 1, . . . , j Ω l ,Υm , the vector n lm,k denotes the outward unit normal vector to S lm,k and the functions g lm,k are defined via the coupling conditions (2.15).

Lumped hydraulic circuits in Υ m . Let the dynamics in each lumped hydraulic circuit Υ m , for m ∈ M = {1, . . . , M }, be described by the vector

y m = [y m,1 , y m,2 , • • • , y m,dm ]
T of state variables satisfying the following system of (possibly nonlinear) ODE system:

dy m dt = A m (y m , t)y m + r m (y m , t) (2.13)
equipped with the initial conditions

y m (t = 0) = y m,0 (2.14)
where y m and r m are d m -dimensional vector-valued functions and A m is a d m × d m tensor. The tensor A m embodies topology and physics of the connections among the circuit nodes and the vector-valued function r m comprises two contributions: (i) sources and sinks within the circuit, including generators of current and voltage; and (ii) connections with Stokes regions.

Coupling conditions. A domain Ω l is connected to a lumped circuit Υ m via the interfaces S lm,k , with k = 1, . . . , j Ω l ,Υm as indicated in Figure 2.9, where we imposed the condition:

g lm,k (x, t) = -P lm,k (t)n lm,k (x) for x ∈ S lm,k and t ∈ (0, T ), (2.15)
where P lm,k is the pressure at the node of the circuit sitting on S lm,k . Under some geometric assumptions on the domain, this condition corresponds to imposing that the average pressure on the interface S lm,k is equal to the nodal pressure, see [START_REF] Heywood | Artificial boundaries and flux and pressure conditions for the incompressible Navier-Stokes equations[END_REF][START_REF] Fouchet-Incaux | Artificial boundaries and formulations for the incompressible Navier-Stokes equations: applications to air and blood flows[END_REF]. In addition, the continuity of mass, and consequently flow rate, across S lm,k implies that

Q lm,k (t) = � S lm,k v l (x, t) • n lm,k (x)dS lm,k for t ∈ (0, T ). (2.16)
For l ∈ L m and k = 1, . . . , j Ω l ,Υm , each term Q lm,k contributes as source/sink for the circuit Υ m ; thus, it is convenient to rewrite r m in (2.13) as

r m (y m , Q lm,k , P lm,k , t) = s m (y m , t) + b m (Q lm,k , P lm,k , t) (2.17)
where s m (y m , t) represents the contribution of sources and sinks within the circuit (generators of current and voltage) and b m (Q lm,k , P lm,k , t) gathers all contributions from the j Υm Stokes-circuit connections.

Remark 2.3.1. Inn order to clearly elucidate the main rationale behind the proposed splitting algorithm, only resistive connections between Stokes regions and lumped circuits are considered in [START_REF] Carichino | Energy-based operator splitting approach for the time discretization of coupled systems of partial and ordinary differential equations for fluid flows: the Stokes case[END_REF]. Resistive connections are among the most commons in blood flow modeling, see for instance [START_REF] Quarteroni | Analysis of a geometrical multiscale model based on the coupling of ODE and PDE for blood flow simulations[END_REF]. In some applications, though, capacitive and inductive elements might be needed. Capacitive connections are used when the fluid pressure in Ω l influences the fluid flow in Υ m and, simultaneously, the fluid pressure in Υ m influence the fluid flow in Ω l , without having actual fluid flow between Ω l and Υ m . Thus, some portions of the boundary of Ω l must be deformable, leading to a fluid-structure interaction problem that goes beyond the scope of this work and might be considered as future research direction. Inductive connections are used when the regime of interest is such that inertial effects become important. Since in the present work we are neglecting inertial effects by adopting the Stokes equations in each Ω l , we also neglect inertial effects in the connections between the Stokes regions and the lumped circuits. We remark that the particular elements allowed in the connections might lead to different initial problems, as pointed out in [START_REF] Moghadam | A modular numerical method for implicit 0D/3D coupling in cardiovascular finite element simulations[END_REF]. We also remark that resistive, inductive and capacitive elements may all be present in the lumped circuit Υ m . Fully coupled problem. The fully coupled problem consists in finding v l , p l , P lm,k , Q lm,k and y m , for l ∈ L, m ∈ L and k = 1, . . . , j Ω l ,Υm satisfying equations (2.7), (2.8) and (2.13), subject to the coupling conditions (2.12), (2.15) and (2.16), the boundary conditions (2.10) and (2.11), and the initial conditions (2.9) and (2.14).

Ω l S lm,k P lm,k Q lm,k Υ m

CHAPTER 2. REDUCED AND MULTISCALE MODELS

Remark 2.3.2. Due to the complex nonlinear structure of the fully coupled problem, to our knowledge, existence of solutions has been proved only for some particular cases. For example, in [START_REF] Quarteroni | Analysis of a geometrical multiscale model based on the coupling of ODE and PDE for blood flow simulations[END_REF], the authors prove local in time existence of a solution, when the connections are made through bridging regions. In [START_REF] Baffico | Multiscale modeling of the respiratory tract[END_REF], the existence of a strong solution for small data is demonstrated, in the case of the coupling between the Navier-Stokes problem and a single resistor.

In [START_REF] Carichino | Energy-based operator splitting approach for the time discretization of coupled systems of partial and ordinary differential equations for fluid flows: the Stokes case[END_REF] we showed that the full coupled Stokes-0d system satisfies an energy identity that embodies the main mechanisms governing the physics of the system. More precisely, the analysis shows that, in the case of lumped circuits having a structure contributing to the overall energy dissipation and in the absence of external forcing terms, the energy of the full coupled problem will decay in time. This is the fundamental physical and mathematical property that we aim at preserving at the discrete level and that provides the main rationale for our splitting scheme described below.

First-order splitting algorithm for the time discretization of the coupled problem and stability analysis

Our numerical strategy is based on the operator splitting technique, see e.g. [Glo03, Chap. II], which leverages a semi-discretization in time to solve sequentially in separate substeps the PDE systems associated with the Stokes regions and the ODE systems associated with the lumped hydraulic circuits. The most important feature of the scheme is that the substeps are designed so that the energy at the semi-discrete level mirrors the behavior of the energy of the full coupled system, thereby providing unconditional stability to the proposed splitting method. We present next the main steps of the algorithm; for more details about its characteristics and the implementation, see [START_REF] Carichino | Energy-based operator splitting approach for the time discretization of coupled systems of partial and ordinary differential equations for fluid flows: the Stokes case[END_REF]. We emphasize that the version of the method detailed below yields, at most, a first-order accuracy in time, since it includes only two substeps; however, the scheme can be generalized to attain second-order accuracy using symmetrization techniques [Glo03, Chap. VI] or suitable time-extrapolations of quantities of particular interest, in the same spirit as [START_REF] Fenandez | Incremental displacement-correction schemes for incompressible fluid-structure interaction[END_REF][START_REF] Fenandez | Fully decoupled time-marching schemes for incompressible fluid/thin-walled structure interaction[END_REF].

First-order splitting algorithm. Let Δt denote a fixed time step, let t n = nΔt and let ϕ n = ϕ(t = t n ) for any general expression ϕ. Let v 0 l = v l,0 for all l ∈ L and y 0 m = y m,0 for all m ∈ M. Then, for any n ≥ 0 solve:

Step 1 For each l ∈ L, m ∈ M and k = 1, . . . , j Ω l ,Υm , given v n l and y n m , find v l and y m such that

∇ • v l = 0 in Ω l × (t n , t n+1 ) (2.18) ρ ∂ v l ∂t = -∇p l + µΔv l + ρf l in Ω l × (t n , t n+1 ) (2.19) dy m dt = b m (Q lm,k , P lm,k , t) in (t n , t n+1 ) (2.20)
with the initial conditions

v l (x, t n ) = v n l (x) in Ω l (2.21) y m (t n ) = y n m (2.22)
and the boundary conditions (2.32)

v l = 0 on Γ l × (t n , t n+1 ) (2.23) � -p l I + µ∇v l � n l = -p l n l on Σ l × (t n , t n+1 ) (2.
v l (x, t n ) = v n+ 1 2 l (x) in Ω l (2.30) y m (t n ) = y
The stability properties of the algorithm are investigated in the next theorem. The analysis was performed on a simplified problem that, however, still preserves the main difficulties associated with the coupling of PDEs and ODEs considered in this work.

Theorem 2.3.3. Consider the coupled Stokes-0d problem described in Section 2.3.1 and assume that (i) the circuit properties are such that, for all m ∈ M, the tensor A m is constant and B m is positive definite, where B m gathers all the contributions from resistive, capacitive and inductive elements in the lumped circuits; (ii) the are no forcing terms, namely F Ω (t) = F Υ (t) = 0 for all t ≥ 0, where F Ω represents the forcing on the system due to body forces and external pressures acting on the Stokes regions and F Υ the forcing on the system due to generators of current and voltage within the lumped circuits, respectively; (iii) the Stokes-circuit connections are of the type described in Section 2.3.1. Under these assumptions, the algorithm given by Steps 1 and 2 is unconditionally stable.

The key ingredient in the proof is the following: the uncoditional stability of the algorithm is a direct consequence of treating the contributions from the j Υm Stokes-circuit connections, represented by b m (Q l,m , P l,m , t), implicitly in Step1 in Eq. (2.20) (for instance, by means of an implicit Euler scheme for the time-discretization). The coupling conditions are not disrupted, and we let the two substeps communicate via the initial conditions. As a consequence, the proposed splitting technique does not introduce uncontrolled artificial terms in the energy, thus ensuring the stability of the algorithm.

Numerical results and discussion

We evaluated the performances of the first order splitting method proposed in Section 2.3.2 on three illustrative examples that differ by: (i) number of Stokes regions and lumped circuits; (ii) number of connections between Stokes regions and lumped circuits; and (iii) type of elements within the lumped circuits. The equivalent electric circuits are represented in Figures 2.10, 2.11 and 2.12, respectively. We derived explicit solutions for the full coupled problem in each example, against which our (as well as other) numerical methods can be tested and we assessed the convergence properties of the method for different choices of time step, showing that the expected first-order convergence in time is actually achieved. The spatial discretization of the Stokes problem is handled via a standard finite-element approach, the temporal discretization of the Stokes and 0d problems is handled via an implicit Euler method and the computational framework relies on the finite element library Freefem [START_REF] Hecht | New development in FreeFem++[END_REF].

p 1 Ω 1 Σ 1 Γ 1 Γ 1 S 11

NUMERICAL SOLVING OF COUPLED DISTRIBUTED-0D PROBLEMS

55

We present the results obtained on Example 1; for Examples 2 and 3, we refer to [START_REF] Carichino | Energy-based operator splitting approach for the time discretization of coupled systems of partial and ordinary differential equations for fluid flows: the Stokes case[END_REF]. They correspond to the case of nonlinear variable resistance R 12 and capacitance C 12 , but we also tested the method in the case of constant R 12 and C 12 against the exact solution and obtained similar performances as in the nonlinear case. Figure 2.13 displays a comparison between the exact solution and different numerical approximations of physical quantities of interest, namely pressure P 11,1 and flow rate Q 11,1 at the Stokes-circuit interface S 11,1 (upper panel) and pressure y 11 and volume y 12 in the 0d circuit (lower panel). We remark that the nonlinearities in R 12 and C 12 have been treated explicitly, which means that R 12 and C 12 are evaluated at the previous time step. A good agreement is obtained for P 11,1 , y 11 and y 12 for Δt = 0.01, and the error decreases as Δt is reduced. The numerical approximation of Q 11,1 captures the periodicity of the solution even for Δt = 0.01, but the peaks are lower than those exhibited by the exact solution. As Δt is the reduced, the computed Q 11,1 approaches the exact peaks of the interface flow rate, thereby capturing the full dynamics of the problem. Figure 2.13 also shows that the numerical solution is not affected by spurious oscillations or instabilities, even for the largest time step. These findings confirm that the choice of the time step affects the accuracy of the computed solution but not the stability of the numerical scheme, as proved in Theorem 2.3.3. Recall that the time-discretization scheme is based on a first order operator splitting technique, which is known to be first-order accurate in time. We performed a standard time refinement study that confirmed this property, as shown in Figure 2 Perspectives. The energy-based operator splitting approach developed in [START_REF] Carichino | Energy-based operator splitting approach for the time discretization of coupled systems of partial and ordinary differential equations for fluid flows: the Stokes case[END_REF] is based on a very simple version of this technique, which is first-order accurate. It has the advantages of an easy implementation, good stability and robustness properties, and low computational cost. This framework could serve as a basis for the extension to high order methods, as briefly discussed in Section 2.3.2. The proposed algorithm allows us to associate to each block of the problem a specific operator and to solve it in separate substeps, by means of a well-suited numerical method. This feature can be further exploited in more complex situations: (i) coupling between Navier-Stokes equations and a lumped parameter model, in the Newtonian or non-Newtonian case; (ii) coupling between a model describing flow in a compliant tube and a 0d appropriate model; (iii) more generally, coupling between different PDEs in the various Ω l , (for instance to accommodate Stokes and Darcy) and reduced models.

Interface quantities

Conclusions and outlook

In this chapter, our strategy has been to tackle the complexity of the bio-fluid dynamics in the coupled eye-cerebral system by means of a two-stage approach. First, we developed network-based models that allow to understand the main dynamics of the system thanks to their simplified mathematical structure and to compute average values of the unknowns at relatively low computational costs. They are also particularly well suited for sensitivity analysis studies to input data, that we successfully used to discriminate the contribution of different factors on flow quantities. In addition, these robust mathematical models have the advantage that can be dynamically complexified and implemented as our understanding evolves. Second, we envisioned a multiscale approach, where regions of interest that require the use of a detailed model are coupled to a reduced 0d description of the rest of the network, thus leading to a coupled PDE-ODE system. We proposed a novel technique based on operator splitting for the time discretization of such coupled distributed-reduced problems, that allows to solve separately and sequentially each model, without the need of sub-iterations.

We selected in the sequel some possible further developments: Towards a multi-domain and multiphysics model of the eye. To the best of our knowledge, the mathematical and computational modeling of the eye is still at its early stages, see for instance [START_REF] Aletti | Mathematical Modelling and Simulations of the Hemodynamics in the eye[END_REF] and the references therein. In the framework of the thesis of Lorenzo Sala, co-advised with Chirstophe Prud'homme (Univ. de Strasbourg) and Giovanna Guidoboni (Univ. of Missouri), we made the first steps in constructing such a full model, in which ocular biomechanics, fluid dynamics and hemodynamics are coupled. Currently, the system combines: (i) a circuit-based model for blood flow in the retinal vasculature, central retinal artery and central retinal vein, as in Sections 2.1 and 2.2.2; (ii) a three-dimensional porous media model for the perfusion of the lamina cribrosa and (iii) a three-dimensional isotropic elastic model for the biomechanics in the lamina cribrosa, retina, choroid, sclera and cornea. From a numerical viewpoint, the problem requires solving a poroelastic model for the coupling between biomechanics and hemodynamics in the lamina cribrosa, which calls for (i) high accuracy in the approximation for both primal variables (i.e. displacement and pressure) and dual variables (i.e. stress and perfusion velocity), and (ii) integral boundary conditions to account for the coupling between zero-and three-dimensional model components. To numerically solve the coupled system, we used a hybridizable discontinuous Galerkin method [START_REF] Cockburn | Unified hybridization of discontinuous Galerkin, mixed, and continuous Galerkin methods for second order elliptic problems[END_REF], implemented in Feel++. Furthermore, to handle the multiscale nature of this problem, the complex coupling has been achieved using an integral boundary condition and a time-splitting energy-based scheme in the spirit of the method developed in Section 2.3. A collaborative effort in this sense was presented and published as peer-reviewed conference abstracts [15; 16] and a full paper is in preparation. Future developments will include the coupling with mathematical and computational models for the cerebro-spinal fluid circulation (see also the discussion in Section 1.6), an extended blood circulation for the whole body and a three-dimensional description of the fluid dynamics of the vitreous and aqueous humor.

Mathematical modeling and simulations for the brain and eye-cerebral system: a multiscale approach. As already discussed in Section 1.6, it would be interesting to extend our contributions on the 3d models for cerebral blood flow to a multiscale framework. Eventually, the 3D description of the flow would be integrated in a model of the entire circulatory system. The advantage of such a closed-loop model are that the artificial boundaries are not present anymore and moreover, it could include feedback and control mechanisms. Similarly, the complexity of the coupled eye-cerebral system calls for a multiscale modeling approach. Network-based models, as the one presented in Section 2.1, allow to capture the main dynamics and averaged behavior, whereas detailed 3d models allow to interface with clinical data that are 3D in nature, e.g. MRA maps or Optical Coherence Tomography (OCT) images. The mathematical complexities of such a multiscale approach, as briefly discussed in Section 2.3, include in particular the extension to 3d geometries and to the case of the Navier-Stokes equations. Also, such an approach requires the development of flexible, versatile and scalable computational framework, in which the multi-components/multiscale/ multiphysics mathematical models can be solved numerically. Moreover, the use of adaptive and robust methods to combine computations with experimental data and measurements is necessary. These new research directions are at the core of the project thesis of Philippe Ricka, co-advised with Christophe Prud'homme (Univ. de Strasbourg) and Giovanna Guidoboni (Univ. of Missouri), starting from September 2017.

Chapter 3

Numerical strategies for coupled fluid-structure problems. Applications to multiphysics computational models in biology and medicine Three-dimensional macro-scale fluid models can describe blood flow dynamics in large vessels (arteries or veins), by means of the Navier-Stokes system, as presented in Chapter 1. Alternatively, in Chapter 2 we used 0d reduced models, based on systems of ODEs describing blood flow with electric analogy, possibly coupled with a 3d description of the flow, in a geometric multiscale approach. The representation of these problems can be enriched to take into account the combined effects of flow and different structures, from a multiphysics perspective. In this direction, two types of fluid-structure interaction phenomena can be considered:

• internal flows, when the external fluid boundary has an elastic part or is entirely composed of an elastic structure, e.g. a model for blood flow in arteries with deformable walls;

• internal structures, when rigid or deformable bodies float or actively move in a fluid flow placed in a container, e.g. a mixture of plasma and red blood cells, as a model for blood at a very small scale.

In this chapter, we gather two contributions, corresponding each to one of these situations. In Section 3.1, a coupled fluid-structure interaction model from the first category is presented, based on a description of blood flow dynamics by the Navier-Stokes equations in a moving domain, coupled with the linear elasticity equations describing the vessel wall deformation. This coupled multiphysics model is subsequently utilized for simulating blood flow in the aorta, in view of a specific biomedical application. The case of an internal structure, namely a rigid particle evolving passively in a fluid, is presented in Section 3.2, with an emphasis on the issue of how to handle contact with the wall of a rigid container.

Computational modeling of blood flow in the aorta

On the basis of [START_REF] Szopos | Computational modeling of blood flow in the aortainsights into eccentric dilatation of the ascending aorta after surgery for coarctation[END_REF], we aim at providing in this section an example of how combining two complementary studies: 60 CHAPTER 3. COUPLED FLUID-STRUCTURE PROBLEMS

• a 3d MRI angiography of the aorta in 30 patients after coarctation repair and

• a 3d fluid-structure computational model characterizing the fluid dynamics in the aorta (flow, pressure, and wall shear stress) can shed some light on the causes of a pathological situation and thus provide interesting insights to the medical community.

Mathematical model and numerical approach

Formulation of the problem. The mathematical model that we consider for blood flow in the aorta is a three-dimensional strongly coupled fluid-structure problem. In the sequel, the superscripts f and s are introduced to identify quantities referring to the fluid and structure problems, respectively. Blood is modeled as a homogeneous, incompressible fluid, with "standard" Newtonian behavior, occupying the moving domain Ω f (t); the dynamics is described by the Navier-Stokes equations in an Arbitrary Lagrangian Eulerian (ALE) formulation, see Equations (3.1-3.3) hereafter. The vessel wall is described as an elastic solid, occupying a domain denoted Ω s (t); the displacement is described in a Lagrangian framework by the linear elasticity equations (3.4-3.6). At the fluidsolid interface Γ w (t), i.e. the common boundary between Ω f (t) and Ω s (t), the coupling conditions (3.7-3.9) are imposed. The boundaries Γ in and Γ out correspond to the artificial inlet and outlet sections, respectively (at the levels where the vessels are truncated). For the geometrical setting and the notations, see also Figure We summarize here the equations for the coupled fluid-structure problem:

Fluid: Coupling conditions:

ρ f � ∂u f ∂t| x 0 + (u f -w) • ∇u f � -2µ∇ • D(u f ) + ∇p = 0, in Ω f (t), (3.1)
∇ • u f = 0, in Ω f (t), (3.2) σ(u f , p)n = g in-out , on Γ in-out . (3.3) Solid: ρ s ∂ 2 d s ∂t 2 -∇ 0 • (F (d s )S(d s )) = 0, in Ω s 0 , (3.4)
d s = 0, on Γ D 0 , (3.5)
geometry: d f = Ext(d s | Γ w 0 ) , in Ω f 0 , Ω f (t) = (Id + d f )(Ω f 0 ), (3.7)
velocity:

u f = w, on Γ w (t), (3.8)
stress:

F (d s )S(d s )n 0 = J(d f )σ(u f , p)F (d f ) -T n 0 , on Γ w 0 . (3.9)
where the unknowns are: the fluid domain displacement d f , the fluid velocity u f and pressure p and the structure displacement d s . For the complete set of notations and more details about the ALE formulation and ALE derivative ∂u f ∂t|x 0

, see [START_REF] Formaggia | of MS&A[END_REF]Chap. 3].

Besides the challenges of three-dimensional realistic simulations for blood flow already discussed in Chapter 1, additional difficulties need to be accounted when dealing with mathematical and computational fluid-structure models per se. As a consequence, a natural question should be considered: is this level of complexity necessary for the problem at stake? In the aorta, the radius may vary in a range of 5 to 10 % between diastole and systole, therefore a realistic model needs to incorporate wall deformation. Moreover, for the clinical application in view, vessel compliance is known to be a key factor for providing more realistic values for the wall shear stress, compared with simulations assuming rigid wall behavior [START_REF] Figueroa | A coupled momentum method for modeling blood flow in three-dimensional deformable arteries[END_REF][START_REF] Ladisa | Aortic coarctation: recent developments in experimental and computational methods to assess treatments for this simple condition[END_REF]. For all these reasons, we adopted a fully coupled fluid-structure model in the sequel.

From the mathematical viewpoint, the analysis of coupled models for fluid-structure interaction raises numerous difficulties because of the nonlinear behavior of the system, which is due not only to the possibly nonlinear equations for each sub-problem, but also to the nonlinearities introduced by the coupling. More information can be found in [START_REF] Formaggia | of MS&A[END_REF]Chap. 8] Discretization. A standard finite-element based method was used for the space discretization. As for the time-advancing scheme, let us briefly discuss the approach used in [START_REF] Szopos | Computational modeling of blood flow in the aortainsights into eccentric dilatation of the ascending aorta after surgery for coarctation[END_REF]. Let Δt be a given time step, at time t n = nΔt, we solve:

(Fluid) Geometry sub-problem: Find the displacement d f,n+1 and an approximation of the fluid domain Ω f (t n+1 ).

Solution: provide d s,n+1 on the interface, then take a harmonic extension.

Fluid sub-problem : Find (u n+1 , p n+1 ) the fluid velocity and the fluid pressure. Solution: implicit Euler scheme.

Solid sub-problem : Find (d s,n+1 , u s,n+1 ) the displacement of the structure and the structure velocity. Solution: mid-point scheme.

Schematically, these steps can be represented as follows:

Geometry + fluid sub-problem: ALE fluid solver

d s,n+1 � → � d f,n+1 , u n+1 , p n+1 � . Solid sub-problem: solid solver � d f,n+1 , u n+1 , p n+1 � � → d s,n+1 .
As a consequence, for each time step, a fixed-point problem needs be solved. At the time when simulations in [START_REF] Szopos | Computational modeling of blood flow in the aortainsights into eccentric dilatation of the ascending aorta after surgery for coarctation[END_REF] were performed, after a thorough investigation of different numerical methods, we chose to solve the nonlinear coupled problem with an implicit algorithm, via Newton's method [FQV09, Chap. 9]. The Jacobian was fully computed [START_REF] Fernández | A Newton method using exact Jacobians for solving fluid-structure coupling[END_REF]. The advantage of such an approach lies in its stability properties (in contrast to explicit schemes, which were proved unstable due to the added mass effect [START_REF] Causin | Added-mass effect in the design of partitioned algorithms for fluid-structure problems[END_REF]), but it leads to an overall high computational cost. Since then, the design of numerical schemes for the fluid-structure problem evolved tremendously: a deep discussion on the topic is beyond the scope of this chapter; we refer to the above-cited reviews and the references therein.

Numerical results and clinical applications

The starting point that motivated [START_REF] Szopos | Computational modeling of blood flow in the aortainsights into eccentric dilatation of the ascending aorta after surgery for coarctation[END_REF] was a question of D. Bonnet and P. Ou (Hôpital Necker) who were interested in a deeper understanding of the underlying bio-mechanical factors involved in the post-operative stages of the coarctation of the aorta. This condition is a congenital defect which consists in a narrowing of the vessel (incidence : 1/100.000 birth), that requires surgical removal of the narrowed segment in the first months of life. The main cause of death for adult patients having undergone a successful repair in infancy is hypertension [START_REF] Clarkson | Results after repair of coarctation of the aorta beyond infancy: a 10 to 28 year follow-up with particular reference to late systemic hypertension[END_REF][START_REF] Simsolo | Long-term systemic hypertension in children after successful repair of coarctation of the aorta[END_REF], but the causes of this condition are still unclear and matter of debate. In a series of previous studies [Ou+04; Ou+06; Ou+07], our clinical collaborators investigated the role of geometry of the aortic arch in resting hypertension through blood pressure measurements at rest and MRI of the aortic arch and left ventricle. Three categories of aortic arch shape were defined: gothic, crenel and normal and the gothic geometry was associated with resting hypertension, as illustrated in Figure 3.2. The limitation in these studies was that despite obvious correlation between the two phenomena, the existence of a causality relationship between abnormalities of aortic arch architecture and high blood pressure at rest could not be established.

From the computational modeling standpoint, several studies [LaD+11a; LaD+11b; Men+] proposed coupled computational and imaging techniques to increase the understanding of coarctationinduced abnormal hemodynamics and vascular biomechanics. However, to the best of our knowledge, at the time of the publication of [START_REF] Szopos | Computational modeling of blood flow in the aortainsights into eccentric dilatation of the ascending aorta after surgery for coarctation[END_REF]), only one publication [START_REF] Olivieri | Hemodynamic Modeling of Surgically Repaired Coarctation of the Aorta[END_REF] explored the use of a computational fluid dynamics model in patients with an abnormal aortic arch shape after coarctation repair, but without taking into account the interaction with the vessel wall.

We successfully performed numerical simulations of the fluid-structure model presented in Section 3.1.1 in the framework of the open-source finite element library LifeV [Lif], using computational meshes of idealized geometries that we generated in the open-source code Modulef [Mod]. The outcome of our numerical results showed abnormal patterns of flow and pressure in the gothic aortic arch compared to the normal Romanesque arch, as illustrated in Figures 3.3 and 3.4, demonstrating that geometrical changes in the aortic arch shape lead to an alteration of the local fluid dynamics in the ascending aorta. In addition, wall shear stress was significantly greater in the gothic arch compared to the "normal" one, affecting, in particular, the anterior and posterior segments of the These new findings on the wall shear stress, an important force that plays a significant role on vascular modifications such as aortic dilatation [Fry+08; Bie+11], but which cannot be directly measured by imaging techniques, were considered of high interest by our colleagues medical doctors. Therefore, they conducted a second complementary 3D-MRI angiography study, that confirmed similar findings, by showing in vivo that the gothic arch is associated with dilatation that is eccentric rather than concentric. Eventually, the results of our two complementary studies were published in a medical journal [START_REF] Szopos | Computational modeling of blood flow in the aortainsights into eccentric dilatation of the ascending aorta after surgery for coarctation[END_REF] and had a direct impact in routine clinical practice: a recommendation for patients with a gothic arch to be seen every 1 to 2 years if increased intimamedia thickening and/or eccentric dilatation of the aorta are present and every 2 years in the absence of such abnormalities; in comparison, patients with a normally smooth Romanesque arch could be evaluated less frequently (every 2 to 3 years). Discussion and outlook. In this computational modeling work, we selected two idealized shapes for the aorta, rather than studying a variety of shapes and angulations. Both geometry and boundary data were not patient-specific, and were not taking full advantage of the MRI complementary study conducted by the clinicians. In addition, a model derived from patient-specific data would have included the anastomotic site, which will usually be less compliant than the normal aorta and, thus, could, itself, contribute to the abnormal flow dynamics. Notwithstanding these limitations, at the time of [START_REF] Szopos | Computational modeling of blood flow in the aortainsights into eccentric dilatation of the ascending aorta after surgery for coarctation[END_REF], this 3D fluid-structure interaction model represented a significant step forward in understanding the complex behavior of blood flow in the aorta, with meaningful clinical implications. Building upon the new theoretical and computational tools developed in the framework of AngioTK [Ang], as presented in Chap. 1, the geometric description, the biomechanical model and the numerical approach could be improved.

From a different perspective, understanding that arch morphology plays an important role in the alteration of the local fluid dynamics lead us to the question of a possible interpretation of the problem in terms of shape optimization. Several works reviewed in [START_REF] Marsden | Optimization in Cardiovascular Modeling[END_REF] have already demonstrated the potential of coupling a shape optimization approach to cardiovascular blood flow simulations with the aim of improving surgical design. In current discussions with Y. Privat (Univ. Paris 6) and his PhD student F. Omnès, we are exploring the possibility of numerically searching the optimal shape of the aorta by minimizing the dissipated viscous energy of the fluid, based on a theoretical and computational framework recently developed in [START_REF] Dapogny | Geometrical shape optimization in fluid mechanics using Freefem++[END_REF]. This research direction will hopefully benefit from the advances of the research project Analysis and simulation of optimal shapesapplication to life sciences, supported by the Paris City Hall (PI: Yannick Privat, dec. 2016-dec. 2020), in which I'm involved as an associated collaborator.

CHAPTER 3. COUPLED FLUID-STRUCTURE PROBLEMS

Discretization in time in simulations of particulate flows

Models and simulations of the coupled dynamics of rigid or deformable bodies evolving in a surrounding fluid are of interest in several biological applications, as for instance understanding the microcirculation of blood [START_REF] Rahimian | Petascale direct numerical simulation of blood flow on 200k cores and heterogeneous architectures[END_REF] or swimming bacteria [START_REF] Decoene | Microscopic modelling of active bacterial suspensions[END_REF]. Among different approaches proposed in the literature, we focus in Section 3.2.1 on a microscopic model, describing the dynamics of rigid particles in a Stokes flow. Our methodological contribution [START_REF] Hillairet | On discretization in time in simulations of particulate flows[END_REF] tackles the issue of handling contacts (either between two particles or between a particle and a wall), which is a key ingredient in the direct numerical simulation of particulate flows, especially in the dense regime.

Dense fluid-particle flows and difficulties arising in their simulations

One of the challenges for fluid/particle flow simulations is to provide an accurate resolution of the so-called lubrication regime, when the distance between two particles immersed in the fluid or between a particle and the wall becomes very small. The presence of extremely high gradients of the velocity in the narrow gap between the particle and the wall (or between the two particles) results in very strong drag forces (which can be referred to as lubrication forces in this case) and thus requires special attention in order to construct efficient and accurate discretization methods.

We consider in the sequel the case of one rigid particle immersed in a viscous, incompressible fluid with a particular emphasis on situations when the particle approaches a plane. Despite its apparent simplicity, this situation encompasses the main difficulties of the lubrication regime and motivated the work in [START_REF] Hillairet | On discretization in time in simulations of particulate flows[END_REF]. The fluid (with the particle inside) fills a fixed domain Ω⊂R d with d = 2 or 3, and the region occupied by the particle B t ⊂ Ω varies with time t, see Figure 3 where D(u) stands for the symmetric gradient of u and n is the unit normal vector on ∂B t directed towards the fluid domain. Using these notations we write out the equations of motion of the particle as follows � m dV dt = F(B t , V, ω) + mg, I t dω dt + ω × I t ω = T(B t , V, ω), (3.12) where m is the mass of the particle and I t is its inertia tensor, expressed in the fixed Cartesian frame and thus dependent on time. Equations (3.12) are coupled with the equations describing the propagation of the particle:

� Ġ = V, ṙi = ω × r i , i = 1, . . . , d. (3.13)
The force F is a sum of the Archimedes force due to gravity and of the drag force which is purely hydrodynamic, i.e. obtained from (3.10) by setting g = 0. The particularity of the drag is that it tends very rapidly to infinity when the particle approaches the wall, thus preventing collisions between them. Indeed, it has been proved in [START_REF] Hillairet | Lack of collision between solid bodies in a 2D incompressible viscous flow[END_REF] (2D case), [START_REF] Hillairet | Collisions in three-dimensional fluid structure interaction problems[END_REF] (3D case) that a smooth rigid body embedded in a viscous fluid cannot touch the wall in finite time. For an overview of existence and uniqueness results for the coupled fluid-structure interaction problem described by Equations (3.10-3.13), including the (near-)contact regime, see [START_REF] Hillairet | Topics in the Mathematical Theory of Interactions of Incompressible Viscous Fluid with Rigid Bodies[END_REF].

In the regime of very small distance between the particle and the wall, the lubrication force is very difficult to take into account in a numerical simulation [START_REF] Lefebvre | Numerical simulation of gluey particles[END_REF]. The corrections specific to the lubrication regime were proposed in the context of the boundary integral methods [START_REF] Qi | Effective moduli of particulate solids: The completed double layer boundary element method with lubrication approximation[END_REF][START_REF] Mammoli | The treatment of lubrication forces in boundary integral equations[END_REF], in that of the Stokesian Dynamics [BB88; SB01], and in that of the force-coupling method [START_REF] Dance | Incorporation of lubrication effects into the forcecoupling method for particulate two-phase flow[END_REF]. Short-range repulsion forces mimicking the lubrication ones were also used in [START_REF] Glowinski | A distributed Lagrange multiplier/fictitious domain method for particulate flows[END_REF] in the context of fictitious domain finite element methods. However, the influence of these (not necessarily realistic) forces on the accuracy of a simulation is not well understood. Another simple idea is just to stop the particle when it tries to penetrate the wall during a numerical simulation, but it is then not necessarily clear which criterion should be chosen to decide if the particle should eventually bounce off the wall and when it should it happen. These questions have a partial answer in the articles [Mau07; Lef09] on the gluey particle model. It is shown there that the particle trajectory satisfies an integro-differential equation in the limit of vanishing viscosity, which is easy to discretize in time using moderate time steps and which predicts the moment of an eventual rebound from the wall. We pursued a similar idea in our approach described in the next section, but our aim is to construct an approximated trajectory of the particle in the lubrication regime that would be accurate enough for any given value of the viscosity, not necessarily small.

A model ordinary differential equation with lubrication forces

Formulation of the problem and construction of an approximated solution. To describe • we first assume that q(t) is the same as q(t) until the latter hits q s for the first time at t = t 1 ;

• next, the trajectory q(t) is frozen until the time t = t2 , which is the first instant after t 1 when v(t) = 0, where v(t) = q(t 1 ) + � t t 1 g(s)ds.

(3.16)

• afterwards, the trajectory resumes again as a solution to (3.15) starting from q s with zero velocity: q(t) =    q(t), for 0 < t < t 1 q s , for t 1 ≤ t < t2 solution to (3.15) with q( t2 ) = q s , q( t2 ) = 0, for t ≥ t2

(3.17)

In the in the special case when g(t) is given by g(t) = � g -(t) < 0, for t ≤ t 0 , g + > 0, for t > t 0 ,

with some positive constants t 0 , g + and a negative function g -(t), with the above notations, we prove the following result: Proposition 3.2.1. If g(t) is given by (3.18), then

0 ≤ t 2 -t2 ≤ 1 n(q s ) , sup t∈[0, t2 ] |q(t) -q(t)| ≤ q s , (3.19)
| q( t2 ) -q( t2 )| = q( t2 ) ≤ g + n(q s ) .

Remark 3.2.3. Proposition 3.2.1 should be interpreted as an estimate of the approximation error and also as a guideline on how to choose the threshold value q s , depending on g and the time step in such a way that the quantity g + n(qs) has the same order of magnitude as the precision in the numerical scheme. Since n(q) → 0 when q → 0, the smaller the threshold, the smaller the error. However, taking q s too small may deteriorate the accuracy of the solution (indeed, the original goal of the introduction of q s was to avoid the extremely small values of q). On the other hand, taking q s too large may unnecessarily perturb the solution in the regions where a straightforward discretization could give more precise results.

Numerical strategy and results. We analyzed the performances of three different algorithms to solve (3.15). Algorithm 1 consists in a straightforward Euler discretization in time of (3.15), but it does not always work in the near-contact regime, since it does not necessarily provide a positive approximation q k . To tackle this issue, we propose in Algorithm 2 the construction of the approximated trajectory q, with an a priori chosen threshold value q s > 0, Finally, we combine in Algorithm 3 an adaptative strategy for the choice of the time step with the threshold approximation given by Algorithm 2.

As an illustration, we show in Fig. 3.7 a comparison between Algorithm 1, unable to produce a physically acceptable solution even with very small Δt = 0.0001 and Algorithm 2, which provides satisfactory results, in a two-dimensional setting for ε = 10 -3 . In order to better observe the quasicontact region and the effect of introducing the threshold in Algorithm 2, we zoom on small distances by passing to a log-scale in Figure 3.8 (left panel). The results are compared with computations using Algorithm 3 with adapted stepsize and taking a number of values for Δt min (right panel). We give in particular the results with Δt min = 0, so that Algorithm 3 is reduced to a standard time marching scheme with automatically adapted step sizes. In practice, we noticed that step sizes Δt of the order 10 -4 were actually sufficient to satisfy our error tolerance criterion almost everywhere apart from a very small range of time around 1, where Algorithm 3 chooses Δt of the order 10 -8 . Setting Δt min to 10 -3 or 10 -4 avoids such small time steps and gives fairly good results. with n(q) = ε/q 3/2 , ε = 10 -3 (log scale). Left: solution obtained by Algorithm 2, choosing q s so that n(q s ) = 1/(20Δt). Right: results obtained by Algorithm 3 with tol = 10 -5 and different values of Δt min specified in the legend.

Towards more complex configurations. Although developments in [START_REF] Hillairet | On discretization in time in simulations of particulate flows[END_REF] are only for the model equation (3.15), we have in mind a possible incorporation of our time marching schemes into a full simulation of flow with solid particles. The coefficient n(q) will be then given by a solution to the Stokes system at each time step where q represents the surface-wall distance in the case of a particle approaching to the wall or the inter-particle distance in the case of two approaching particles. The "external" force g(t) in (3.15) will then include the components of the physical external forces, like gravitational force, in the direction connecting the closest points on the surfaces of nearly touching particles, plus other corrections. Most importantly, these include the part of hydrodynamic drag due to the convective terms in the Navier-Stokes system, if the convective terms in the fluid motion equation are not neglected. The time marching scheme should thus work robustly for any right hand side g(t), in particular for a function changing the sign. Indeed, a possible scenario even in the simplest test case of a particle falling under its weight on a horizontal wall (so that g(t) is negative in the beginning) could be that the convective terms in Navier-Stokes equations become so important at some point in time that they prevail over the gravitation force and g(t) becomes positive, which may lead to the rebound of the particle. The next step could take into account the case of deformable bodies evolving in a fluid flow.

Another perspective of the present approach is to include more complex solid movements. In this direction, it would be interesting to carry out simulations in the following configurations: (i) kissing phase in the "Drafting-Kissing-Tumbling" [START_REF] Fortes | Nonlinear mechanics of fluidization of beds of spherical particles[END_REF]: the fall of a homogeneous disk under the action of the gravity force in a cavity which contains a hole in form of a disk; (ii) the sedimentation of a particle which is almost in contact with a vertical wall [START_REF] Cox | The motion of suspended particles almost in contact[END_REF]. In both cases the problem would be solved in the lubrication approximation.

Conclusions and outlook

We summarized in this chapter two contributions where biomedical questions were treated from the mathematical and computational viewpoint as coupled fluid-structure interaction problems. Based on a careful description of the bio-mechanical phenomena into play, we addressed these questions by means of efficient numerical methods, either selected within the state-of-the-art context or newly developed for the specific problem in view. In the previous chapters, we focused on the fluid mechanical aspects of the problems, governed by the Navier-Stokes equations or by reduced fluid models, possibly coupled in a geometrical multiscale framework. Here the complexity and the challenges came from the multiphysics approach, in which the fluid equations are coupled with a system of equations describing the dynamics of a rigid or deformable solid.

Chronologically, these were the first works I developed in mathematical biology and numerous perspectives are now possible due to recent developments in the field. Besides the extensions already mentioned at the end of each section, a last comment, taking as a starting point the examples treated in this chapter. On one hand, the two problems share a similar mathematical structure and therefore a unified language was utilized to solve them, with new analytical and numerical methods. On the other hand, the first one involves a complex three-dimensional, highly nonlinear fluid-structure model, with the goal of answering a medical question, whereas the second one is a toy model that brings a methodological answer to a simplified, yet meaningful biological problem. The two of them are illustrative and somehow dual aspects in the interplay between what we know and what we want to solve, both in mathematics and in biology or medicine, with infinite possibilities of cross-fertilization of ideas. In this spirit, the conclusion of the chapter, and more generally of the manuscript, is the following call for new challenges and opportunities:

Mathematics is biology's next microscope, only better; biology is mathematics' next physics, only better.

Joel Cohen (Rockefeller and Columbia Universities, NY, 2004).

This manuscript gathers my contributions focused on developing new mathematical and computational methods for analyzing biological flows as complex multiphysics and multiscale phenomena. The description of the underlying mechanisms stems from the basic principles of fluid dynamics and is translated into systems of partial or ordinary differential equations. The overall objective of this work is the study of these equations at the continuous and discrete levels, their coupling and the development of a reliable and efficient computational framework to implement various numerical methods to approximate the solutions to these problems. The numerical simulations incorporate realistic geometries, are thoroughly validated against experimental data and target specific biomedical applications. The first chapter focuses on three-dimensional models, in which the motion of a biofluid in a complex, realistic geometry is governed by the incompressible Navier-Stokes equations in a domain with inlet and outlet boundaries, the main application in view being the study of the cerebral venous network

Introduction 1

 1 Three-dimensional models, analysis and efficient simulations for fluid equations. Applications to computational hemodynamics in the cerebral network 1.1 Mathematical modeling of blood flow in the cerebral venous system . . . . . . . . . . 1.1.1 A new three-dimensional model for large cerebral veins . . . . . . . . . . . . . 1.1.2 Numerical results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 A novel formulation of the Stokes system involving pressure boundary conditions . . 1.2.1 General framework and the Lagrange multiplier formulation . . . . . . . . . . 1.2.2 Numerical results and applications . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Influence of different modeling assumptions . . . . . . . . . . . . . . . . . . . . . . . 1.3.1 Sensitivity analysis framework . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Numerical results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Large-scale blood flow simulations and validation . . . . . . . . . . . . . . . . . . . . 1.4.1 Computational framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.2 The FDA benchmark nozzle model . . . . . . . . . . . . . . . . . . . . . . . . 1.5 An open-source framework to generate virtual MRA images from real MRA images . 1.6 Conclusions and outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Reduced and multiscale mathematical and computational models for biofluids.

Figure 1 .

 1 Figure 1.1: 1: internal jugular veins, 2: vein of Galen, 3: straight sinus, 4: confluence of sinuses, 5: lateral sinus (transverse portion), 6: lateral sinus (sigmoid portion), 7: superior sagittal sinus, 8: internal cerebral vein, 9: basilar vein, 10: superior cerebral veins, 11: superior anastomotic veins.

Figure 1 . 2 :

 12 Figure 1.2: Top: cerebral venous network, visualization of the flow patterns: instantaneous streamlines, colored with velocity magnitude. Bottom: zoom on the velocity field in the confluence of sinuses (left panel) and in the superior sagittal sinus (right panel).

Figure 1 . 3 :

 13 Figure 1.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, St. S: Straight Sinus, L. TS: Left Transverse Sinus, R. TS: Right Transverse Sinus.

  and such that each connected component of Γ 2 is flat, represents a partition without overlap of the boundary of Ω and n indicates the outward normal to ∂Ω. The notations are similar to Section 1.1, the function f is a given external force, the function p 0 a given pressure, and d = 2, 3 is such that Ω ⊂ R d .

Figure 1 . 4 :

 14 Figure 1.4: Cerebral venous hemodynamics obtained by imposing a pressure drop between the inlet and outlet sections, pressure field (left panel) and streamlines (right panel).

Figure 1 . 5 :

 15 Figure 1.5: Zoom on velocity vectors at some inlet sections (left panel), respectively outlet sections (right panel).

2 . 1 ;Figure 1 . 6 :

 2116 Figure 1.6: Representation of the relative importance of various factors in models of the different hemodynamic scales: ++ indicates primary importance; -indicates secondary or negligible importance; ?? indicates potential or unclear importance (taken from Creative Commons and [Ste12]).

Figure 1 . 7 :

 17 Figure 1.7: Effects of inlet velocity magnitude v in . The reference solution for each level mesh is computed using v in = 10 -3 m s . Continuous lines correspond to v in = 20 -3 m s and dashed lines to v in = 30 -3 m s .

Figure 1 . 8 :

 18 Figure 1.8: Effects of inlet velocity magnitude v in . The reference solution for each model is computed by using v in = 10 -3 m s . Continuous lines correspond to v in = 20 -3 m s and dashed lines to v in = 30 -3 m s .

Figure 1 . 10 :

 110 Figure 1.10: Re T hroat = 500, experimental data vs. numerical results: normalized axial velocity (left); normalized pressure difference along z (right).

Figure 1 . 11 :

 111 Figure 1.11: Re T hroat = 500: validation metrics E z (left) and E Q (right).

Figure 1 .

 1 Figure 1.12: VIVABRAIN pipeline.

Figure 2 . 1 :

 21 Figure 2.1: Fluid dynamic connections between the eye and the brain: (a) The blood supply to the eye is secured by the ophthalmic artery branching from the internal carotid artery [Eyeb]; (b) The blood drains from the eye into the cavernous sinus, which also receives blood from superficial cortical veins [Eyec]; (c) The cerebrospinal fluid (CSF) fills the subarachnoid space which surrounds the brain and the optic nerve [Eyea].

2 :

 2 Figure 2.2: Network model of fluid flows in the brain and eyes. The nodes correspond to the connection between the brain and eye models. The connection Intracranial Arteries-Ophthalmic Artery represents arterial supply; the connection Venous Sinuses-Cavernous Sinus represents the venous drainage; the grey and light blue arrows represent the pressures acting on both sides of the lamina cribrosa.

Figure 2 . 3 :

 23 Figure 2.3: Network model for the eye vasculature (black portion) and aqueous humor production and drainage (blue portion). The vasculature comprises circulation of blood in retinal, ciliary and choroidal vascular bed. Variable resistances are marked with arrows.

Figure 2 . 4 :

 24 Figure 2.4: Simulated IOP, ICP and laminar compressive stress in the LHDT, M0, M1 and M2 conditions. Dotted lines indicate representative Δπ values in the four simulated cases, which differ from the physiological value on earth of 25mmHg due to fluid shift. Dashed lines indicate Δπ values for which veins collapse in retinal (black), choroidal (red) and ciliary (blue) circulations.

Figure 2 . 5 :

 25 Figure 2.5: Simulated normalized blood flow in retinal, choroid and ciliary circulations. LHDT and M2 cases are selected as most representative. Dotted lines indicate representative Δπ values in the four simulated cases, which differ from the physiological value on earth of 25mmHg due to fluid shift. Dashed lines indicate Δπ values for which veins collapse in retinal (black), choroidal (red) and ciliary (blue) circulations.

Figure 2 . 6 :

 26 Figure 2.6: Probability density function of IOP (left panel) and Sobol indexes (right panel), resulting from the sensitivity analysis performed on the mathematical model of equation (2.1) when variations in ciliary capillary blood pressure cBP , total inflow facility L, blood/aqueous humor osmotic pressure difference Δπ s , trabecular outflow facility C 0 , uveoscleral outflow facility k 1 and episcleral venous pressure EV P are considered.

Figure 2 . 7 :

 27 Figure2.7: Equivalent electrical circuit for aqueous humor production and drainage. Each node of the circuit (black bullet) is connected to another node either through a conductance (denoted by the symbol G) or through a current source (denoted by a rhombus or by a circle). The nodal quantities indicated by P (capital letter) represent a pressure (in mmHg). The nodal quantities indicated by C (capital letter) represent a concentration (in g ml -1 ) while those indicated by m (small letter) represent a molar concentration (in µmol ml -1 ). The quantities surrounded by a box are assumed to be given values defined by the user. The current sources I 1 and I 2 are fixed values while the other current sources are controlled by the drop of the two nodal values across which they are located or by the drop of the two nodal values to which they are connected by the dashed arrows. The vertical shaded triangle symbol represents connection of the considered network branch to ground.

Figure 2 . 8 :

 28 Figure 2.8: Sensitivity analysis as a function of percentage change in D AH . Stromal and anterior chamber pressures are influenced by variations in aqueous humor drainage, while intraocular pressure is influenced in a milder manner. Each pressure tends to a saturation value as drainage is increased. A red square indicates the baseline value for each pressure.

Figure 2 . 9 :

 29 Figure 2.9: Schematic representation of the coupling between the Stokes region Ω l and the lumped circuit Υ m . Coupling conditions for the pressure P lm,k and the flow rate Q lm,k should be imposed on the interface S lm,k .

=

  v l (x, t n+1 ) and y n+1 m = y m (t n+1 ).

Figure 2 . 13 :

 213 Figure 2.13: Example 1.Comparison between the exact solution and the corresponding numerical approximation for interface quantities (top) and 0d unknowns (bottom), for three time steps Δt = 0.01, 0.005, 0.001, over one period once the periodicity is reached.

Figure 2 . 14 :

 214 Figure 2.14: Plot of the energy norm errors, in logarithmic scale as a function of the global time step Δt = 0.01, 0.005, 0.001 for the three examples considered.

  3.1.

Figure 3 . 1 :

 31 Figure 3.1: Geometrical setting and notations for the coupled fluid-structure problem.

3. 1 .F

 1 COMPUTATIONAL MODELING OF BLOOD FLOW IN THE AORTA 61 (d s )S(d s )n s 0 = 0 on Γ N 0 . (3.6)

  or [BGN14, Chap. 1 and Chap. 2]. The numerical treatment of the problem of course inherits these difficulties and a lot of research work has been devoted to propose efficient algorithms to solve it, as reviewed for instance in [FQV09, Chap. 9], [BGN14, Chap. 1 and Chap. 2] or very recently in [QMV17].

Figure 3 . 2 :

 32 Figure 3.2: Aortic arch geometry and associated hypertension (HT), from [Ou+06].

Figure 3 . 3 :

 33 Figure 3.3: Flow rate evolution in time in the ascending aorta (red line) and descending aorta (green line); normal arch (left panel), gothic arch (right panel). Note the reverse flow in the gothic arch in late systole.

Figure 3 . 4 :

 34 Figure 3.4: Pressure wave propagation. Top: snapshots corresponding to the moment when the pressure wave reaches the descending aorta, time ≈ 0.016 sec. in the normal and gothic arch. Colors denote pressure values, from smallest (blue) to largest (yellow). Cross-section 1, pressure magnitude plotted against time at 4 points located at the anterior, posterior, and left and right wall in a cross section of the ascending aorta. Cross-section 2, pressure magnitude plotted against time at 4 points located at the superior, inferior, and left and right wall in a cross section of the horizontal aorta (left panels: Romanesque arch; right panels: gothic arch).

Figure 3 . 5 :

 35 Figure 3.5: Wall shear stress (WSS) mapping and magnitude. Top: WSS distribution at an instant corresponding to the passage of the pressure wave (time ≈ 0.007 second). Colors denote WSS values, from smallest (blue) to largest (red). Cross-section 1, magnitude of the WSS plotted against time at 4 points located at the anterior, posterior, and left and right wall in a cross section of the ascending aorta. Cross section 2, magnitude of the WSS plotted against time at 4 points located at the superior, inferior, and left and right wall in a cross section of the horizontal aorta (left panels: Romanesque arch; right panels: gothic arch).

  .

6 .

 6 Assuming that the inertial effects are negligible in the fluid and the no-slip conditions are valid on the boundaries of Ω and B t , the fluid motion is governed by the Stokes equations:

Figure 3 . 6 :

 36 Figure 3.6: Geometrical setting and notations.

Figure 3 . 7 :

 37 Figure 3.7: Solution to the ODE (3.15) with n(q) = ε/q 3/2 , ε = 10 -3 . Left: the exact solution and solutions obtained by Algorithm 1 with Δt ranging from 0.1 to 0.0001. Right: solutions obtained by Algorithm 2 with Δt in the same range. The results obtained with Δt ≤ 0.01 on the right are visually indistinguishable from the exact solution.

Figure 3 . 8 :

 38 Figure 3.8: Solution to the ODE(3.15) with n(q) = ε/q 3/2 , ε = 10 -3 (log scale). Left: solution obtained by Algorithm 2, choosing q s so that n(q s ) = 1/(20Δt). Right: results obtained by Algorithm 3 with tol = 10 -5 and different values of Δt min specified in the legend.

  . The purpose of Chapter 2 is two-fold: (i) first, we present contributions towards the elaboration of several reduced 0d models describing the coupled dynamics of different biofluids in the eyecerebral system; (ii) second, we describe a new splitting strategy for the numerical solving of coupled systems of partial and ordinary differential equations for fluid flows. In Chapter 3, the fluid dynamics description from the previous chapters is enriched to take into account the combined effects of flow and different structures, from a multiphysics perspective. The resulting framework is subsequently utilized for simulating blood flow in the aorta, in view of a specific biomedical application and for the numerical simulation of particulate flows, with an emphasis on the issue of how to handle contacts. INSTITUT DE RECHERCHE MATHÉMATIQUE AVANCÉE UMR 7501 Université de Strasbourg et CNRS 7 Rue René Descartes 67 084 STRASBOURG CEDEX Tél. 03 68 85 01 29 Fax 03 68 85 03 28 www-irma.u-strasbg.fr irma@math.unistra.fr IRMA 2017/007 http ://tel.archives-ouvertes.fr/tel-01646867 ISSN 0755-3390 Ins titut de R echerche Mathématique Avancée

  2.1.1 Construction of the model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.2 Microgravity conditions modeling and numerical results . . . . . . . . . . . . 2.2 A zoom on mathematical models for the aqueous humor . . . . . . . . . . . . . . . . 2.2.1 Modeling the intraocular pressure under uncertainty . . . . . . . . . . . . . . 2.2.2 Modeling the coupled dynamics of ocular blood flow and aqueous humor . . . 2.3 A new algorithm for the numerical solving of coupled distributed-0d problems . . . 2.3.1 Description of the coupled Stokes-0d problem . . . . . . . . . . . . . . . . . . 2.3.2 First-order splitting algorithm for the time discretization of the coupled problem and stability analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.3 Numerical results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Conclusions and outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Numerical strategies for coupled fluid-structure problems. Applications to multiphysics computational models in biology and medicine 3.1 Computational modeling of blood flow in the aorta . . . . . . . . . . . . . . . . . . . 3.1.1 Mathematical model and numerical approach . . . . . . . . . . . . . . . . . . 3.1.2 Numerical results and clinical applications . . . . . . . . . . . . . . . . . . . . 3.2 Discretization in time in simulations of particulate flows . . . . . . . . . . . . . . . . 3.2.1 Dense fluid-particle flows and difficulties arising in their simulations . . . . . 3.2.2 A model ordinary differential equation with lubrication forces . . . . . . . . . 3.3 Conclusions and outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	ii	CONTENTS
	3 References of the author
	Other references
	i

  We carefully performed a dimensional analysis of the model and computed values for the Reynolds, Stokes, and Strouhal number, reported in Table1.1. 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 analysis shows that flow in the cerebral veins is governed by the Navier-Stokes equations for an incompressible viscous fluid, in a quasi-steady regime:

		CHAPTER 1. THREE DIMENSIONAL MODELS
	an appropriate biomechanical model.
	Parameter	Value	Unit	Source
	Mass density ρ Dynamic viscosity µ Cross-sectional velocity V * (clinical studies, jugular vein) Cross-sectional velocity V * (clinical studies, superior sagittal sinus) Reynolds number (approx. values)	1055 3.5 •10 -3 (10 -11) • 10 -2 (8.5 -11.3) • 10 -2 (30 -50) • 10 -2 15 • 10 -2 (15.2 ± 3) • 10 -2	[ kg m 3 ] [P a • s] [ m s ] [ m s ] [ m s ] [ m s ] [ m s ]	[Thi08] [Thi08] [HMH13] [Ogo+11] [Sch04] [Sch04] [Gid+96] [11]
	Internal jugular vein Superior sagittal sinus Stokes number (approx. values)	90 -232 90 -144	--	[11]
	Internal jugular vein Superior sagittal sinus Strouhal number (approx. values)	1.38 -3.84 1.10 -1.75	--	[11]
	Internal jugular vein Superior sagittal sinus	0.014 -0.030 0.013 -0.021	--
	Table 1.1: Values of the mechanical and dimensionless parameters in the model.

Table 1 .

 1 2: Mesh convergence tests; h: characteristic element size, N elt : number of tetrahedra, N dof : number of degrees of freedom.

		h	N elt	N dof
	m0	0.03	322 013 1 640 236
	m1	0.02	985 484 4 717 123
	m2 0.015 2 008 757 9 171 904

Table 1 .

 1 . Strategy P GASM couples monolithically a Krylov iterative solver with an additive Schwarz preconditioner, whereas Strategies P Block Krylov iterative solver with a block preconditioning strategy, following[START_REF] Elman | Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics[END_REF]. To apply this framework, we remark that we deal with a double saddle point problem: we gather either the velocity-pressure unknowns or the velocity-Lagrange multiplier unknowns to setup a two level preconditioner. The comparison clearly shows the limits of Strategy P GASM and the gain in terms of computational time when choosing Strategy P Block

	1	and P Block 2

3: Time comparison for three preconditioning strategies (in seconds). In brackets, the number of iteration used by solver.

Simulations were performed with 96 processors and the time is measured in seconds.

Table 2 .

 2 1: Analogy between electric and hydraulic networks, adapted from [FQV09, Chap.10].

	34	CHAPTER 2. REDUCED AND MULTISCALE MODELS
		Hydraulic network	Electric network
		Pressure	Voltage
		Flow rate	Current
		Volume	Charge
		Blood viscosity	Resistance R
		Blood inertia	Inductance L
		Wall compliance	Capacitance C

  C aq , P cil c , P st , P ant , P post , m bl , P

	O2 c , P O2 v , P O2 ce , P O2 aq	� T	,	(2.2)

  l (x, t n ) • n lm,k (x, t n ) dS lm,k = Q lm,k (t) in (t n , t n+1 ) For each l ∈ L, m ∈ M and k = 1, . . . , j Ω l ,Υm , given v

	2.3. NUMERICAL SOLVING OF COUPLED DISTRIBUTED-0D PROBLEMS	53
	with			�
								(2.26)
				S lm,k
	and then set
	v	n+ 1 2 l	= v l (x, t n+1 ), p n+1 l	= p l (x, t n+1 ) and y	n+ 1 2 m	= y m (t n+1 ).	(2.27)
	Step 2 n+ 1 2 l that	and y m , find v l and y m such n+ 1 2
		ρ	∂ v l ∂t	= 0			in Ω l × (t n , t n+1 )	(2.28)
			dy m dt	= A m (y m , t) y m + s m (y m , t)	in (t n , t n+1 )	(2.29)
	with the initial conditions
								24)
	�	-p l I + µ∇v l	�	n lm,k = -P lm,k n l,m	on S lm,k × (t n , t n+1 )	(2.25)

v

  Figure 2.10: Example 1. The two-dimensional Stokes region Ω 1 is connected to the lumped circuit Υ 1 via a resistive element with resistance R 11 .Figure 2.11: Example 2. The two-dimensional Stokes regions Ω 1 and Ω 2 are connected to the lumped circuit Υ 1 via resistive elements with resistance R 11 and R 14 , respectively.Figure 2.12: Example 3. The two-dimensional Stokes region Ω 1 is connected to the lumped circuit Υ 1 via two resistive elements with resistance R 11,1 and R 11,2 , respectively.

													Υ 1
						,1
					P 11,1			y 11				y 12	� p 1
						Q 11,1	R 11 C 11		C 12 (y 1 , t) R 12 (y 1 , t)	R 13
								Υ 1
												p
	Σ 1	Γ 1	S 11,1						y 12		R 13	S 21,1	Γ 2	Σ 2
	p 1	Ω 1 Γ 1	P 11,1 Q 11,1	R 11 C 11	y 11	R 12	L 12 C 13	y 13	R 14	Q 21,1 P 21,1	Γ 2 Ω 2	p 2
		p2											p1
			R 14		S 11,2		Γ 1		S 11,1		R 12	Υ 1
		y 12			P 11,2		Ω 1		P 11,1		y 11
			R 11,2 C 12	Q 11,2		Γ 1	y 13	Q 11,1	R 11,1 C 11
						R 13	L 13

  .14.

								CHAPTER 2. REDUCED AND MULTISCALE MODELS
			��� � �������������					��� � �������������						��� � �������������
			��� �	�������������					��� �	�������������						��� �	�������������
			��� � �������������					��� � �������������						��� � �������������
	��������������������������	�� �� �� �� �� �� �� �� �� �� �� �	��������������	�� ��	��������������������������	�� �� �� �� �� �� �� �� �	��	��������������	��	��	��������������������������	�� �� �� �� �� �� �� �� �� �� � ��	��	��������������	��	��
			(a) Example 1					(b) Example 2						(c) Example 3

Remerciements -Acknowledgments

CHAPTER 3. COUPLED FLUID-STRUCTURE PROBLEMS

our strategy, we consider the following simplified setting: Ω ⊂ R 2 is the half-plane {(x, y), y > 0} and the particle is a disk of radius R. Let moreover g(t) = g(t)e 2 and assume the particle is at rest at the initial time. The x-component of the particle velocity and its angular velocity will then vanish at all time. The position of the particle is fully determined by its distance q from the bottom, as in Figure 3.6. The net force F is the sum of the drag, which is a function of q and V, linear in V, and of the Archimedes force:

where n(q) is the drag coefficient computed by the Stokes equations (3.10). Denoting the ycomponent of the velocity by v, we are thus led to the following differential equations

In [START_REF] Hillairet | On discretization in time in simulations of particulate flows[END_REF] we provide the following derivation of the asymptotic expansion for the lubrication force acting on a disk immersed in a Newtonian fluid and approaching the wall:

Theorem 3.2.1. The viscous drag force satisfies F = -n(q)V where

with ε(q) → 0 when q → 0.

The proof of Theorem 3.2.1 is divided into three steps. First, we use the variational formulation for the solution (u, p) to (3.10) to compute the associated drag coefficient n(q) in the expression for viscous drag F = -n(q)V. Second, we deduce from the variational formulation the lower and upper bounds for n(q). We conclude the proof by computing an asymptotic expansion of the bounds of this range. Remark 3.2.2. The equivalent result in the three-dimensional setting is well-known but to the best of our knowledge it was new in two dimensions. The proof method also extends to the threedimensional setting with small changes and can be adapted to other boundary conditions.

After eliminating v from the system (3.14) and going to non-dimensional variables, we obtain the following equation for q(t): q = -n(q) q + g,

where n(q) = ε/q 3 2 in the two-dimensional case, n(q) = ε/q in the three-dimensional case, respectively and ε depends on the characteristic values of g, of the fluid and the solid. We remind that equations of the type (3.15) are at the basis of the gluey particle model of [START_REF] Maury | A gluey particle model[END_REF][START_REF] Lefebvre | Numerical simulation of gluey particles[END_REF], that consists in considering the limit ε → 0, which is physically the limit of vanishing viscosity. Note however that in practice ε computed from the non-dimensionalization procedure is not necessarily small.

In the lubrication regime of q → 0, the drag coefficient n(q) tends to infinity, which can create difficulties when computing numerical solutions of (3.15). In order to overcome this problem, one should either use a very small time step when discretizing the equation, but this can become computationally expensive, or propose a strategy to prevent q from becoming too small. The approach we proposed in [START_REF] Hillairet | On discretization in time in simulations of particulate flows[END_REF], inspired by the latter idea, consists in introducing a threshold value q s for the distance q, below which, instead of solving (3.15), we introduce an alternative model and construct a sufficiently accurate "approximate" solution q instead of the "true" solution q. The construction of the approximated trajectory is the following: