Raoul Van Loon

Thomas Podgorski

férents cas, d'abord avec une ou sans valve, pulsation du fluide puis des parois mobiles. Ce type d'écoulement est similaire à celui des globules rouges dans un flux sanguin.

English -This thesis aims at numerically modelling the lymphatic network at different scales. A failure of this network results in an accumulation of fluid in the area concerned, which can lead to lymphoedema.

With the objective to understand this network and how it develops its capacity to move lymph, we are interested in a comparison between different description and models on the lymphatic network from the literature. We first study the discrete approach in zero dimension and in one dimension. In addition, the different constitutive equations observed in the literature are detailed and analysed. In order to understand the interactions between the basic motile elements of the lymphatic network, different methods for coupling calculations in fluid-structure simulation are presented. Then, with the help of different recent articles, we compare different approaches and geometries for the study of a lymphangion.

In this zero-dimensional approach, a new numerical formulation is used for the calculation of lymph flow in the collecting network. Equations of this model are detailed in this document, and the number of parameters generally used in the constitutive equation is reduced. Moreover, these equations allow for a variable contraction frequency depending on the load imposed by the boundary condition applied. Different specific cases such as divergent and convergent bifurcations, elementary units of a network are first studied. The results from the simulations are compared with experimental data. Finally, a specific and realistic network geometry extracted from an anatomical drawing of a leg is used to simulate the model and show complex synchronization behaviours between lymphangions.

For the last chapter, a two-dimensional model of a lymphangion is proposed, which will then be used to study the behaviour of valves, lymphocytes and walls. The operation of the multi-physics fluid-structure code is explained, it is based on a method called: Immersed Structural Potential Method (ISPM). Initially, the equations of fluid and solid mechanics are introduced, how they are coupled, as well as the details of their implementation. Then, we study the behaviour of a group of lymphocytes in the lymphatic channel using a fluidstructure interaction code. Using the geometry of a lymphangion, we compare the displacement of the lymphocytes in different cases, first with or without valve, pulsation of the fluid and then the moving walls. This type of flow is similar to that of red blood cells in a blood stream.

Summary

Français -Cette thèse porte sur la modélisation numérique du réseau lymphatique à différentes échelles. Une défaillance de ce réseau résulte en une accumulation du fluide dans la zone concernée pouvant entraîner un lymphoedème.

Afin de comprendre ce réseau et comment il développe sa capacité à déplacer la lymphe, nous nous s'intéressons dans un premier temps à une comparaison des différentes description et modèles de réseau lymphatique dans la littérature. Nous étudions d'abord les approches discrètes à zéro dimension et à une dimension. De plus, les différentes équations constitutives observées dans la littérature sont extraites et analysées. Dans le but de comprendre les interactions entre les éléments motiles de base du réseau lymphatique, nous présentons différentes méthodes pour les calculs de couplage en simulation fluide-structure. Puis, en nous appuyant sur différents articles récents, nous comparons différentes approches et géométries sur l'étude des lymphangions (unités fonctionnelles valvulées des réseaux lymphatiques).

Pour cette approche en zéro dimension, une nouvelle formulation numérique est employée pour le calcul d'écoulement de lymphe dans le réseau collecteur. La formulation analytique est détaillée et justifié dans ce document, nous réduisons le nombre de paramètres dans l'équation constitutive généralement utilisé. De plus cette formulation permet d'obtenir une fréquence de contraction variable avec la charge imposée par les conditions limites du système. L'étude de différents cas spécifiques tels que les bifurcations convergentes et divergentes, éléments fondamentaux d'un réseau est effectué. Nous comparons les résultats numériques avec des données expérimentales. Enfin, les résultats d'une géométrie de réseau spécifique basée sur des planches d'anatomie de la jambe sont présentés, montrant des phénomènes de synchronisation complexe entre lymphangions.

Dans le dernier chapitre, un modèle en deux dimensions d'un lymphangion est proposé, qui sera ensuite utilisé pour l'étude du comportement des valves, lymphocytes et parois. Nous expliquons en détail le fonctionnement du code numérique multiphysique d'interaction fluide-structure utilisant une méthode nommée : Immersed Structural Potential Method (ISPM). Initialement, nous introduisons les équations de la mécanique des fluides et des solides, comment elles sont couplées, ainsi que le détail de leurs implémentations. Ensuite, le comportement d'un groupe de lymphocytes dans le canal lymphatique à l'aide d'un code d'interaction fluide-structure est étudié. En utilisant la géométrie d'un lymphangion, nous comparons le déplacement des lymphocytes dans difiii Chapter 1

Introduction

The lymphatic system constitutes a transport system, which was first briefly mentioned by Hippocrates in the fifth century BC. And, first described scientifically only in the seventeenth century by Thomas [START_REF] Bartholin | De lacteis thoracicis in homine brutisque nuperrime observatis, historia anatomica. Publice proposita, Respondente M. Michaele Lysero[END_REF]). The system operates in parallel with the blood system. However by comparison, it does not form a closed loop system, like the blood system with veins, arteries and the heart as a pump. From a functional point of view, it transports excess interstitial fluid back into the blood circulation, via the thoracic duct. It also helps carry some fat from the digestive tract into the blood, or to the appropriate organs. Along with the excess interstitial fluid, some proteins are transported back to the circulation via for example the thoracic duct. The lymphatic system also serves a great purpose for immune cells and contributes to the immune response by carrying lymphocytes for example. In fact, lymph nodes across the network filter the interstitial fluid and break down bacteria, viruses and waste. They also collect antigens carried by the lymph, triggering the release of antibodies and lymphocytes. Furthermore as it moves fluid and cells back into venous system, it contributes importantly to the spread of cancer cells throughout the entire body and tend to cluster in the lymph nodes.

Lymph drainage can be impaired as a consequence of infection, surgery, transplantation, or congenital diseases. This condition is called lymphoedema, in which lymph does not flow properly and tends to swells until it solidifies. As of today, there are no treatments and management for this disease has limited success. As a matter of fact, multiple paths have been followed to limit its effects, showing the limited understanding of initial causes. For example, arm lymphoedema can occur in breast cancer patients. As mentioned earlier, cancer cells tend to cluster in lymph nodes and are often removed which modifies the equilibrium of the lymphatic network and thus its ability to pump.

Despite the importance of the lymphatic system in health and disease, its transport properties have received a rather limited attention. The lymphatic system is complex with different functions between anatomical sites. In fact, mammals have similar lymphatic system, it is present in reptiles and bird as well but with less or no lymph nodes. However, insects have open circulatory fluid, and blood and lymph functions are ensured by a similar fluid called hemolymph.

In this introduction, an overview will be given of the lymphatic system, its structure, function as well as the main pathological issues and clinical implications. The disease impact will be discussed alongside the links to clinical applications. Finally, the outline of this manuscript is described.

The lymphatic system

The lymphatic system, which is used also in the immune system, carries proteins, bacteria and cellular debris. This system is made out of the lymph (lympha: water in Latin), the conducting network and the lymphatics organs. The conducting network connects a fine network (capillary-like) in the tissue and the organs to the lymphatics organs and then to the lymphatic duct. For this study, we will only focus on the lymphatic network, as represented in figure 1.2 which represents the lymphatic network in leg skin.

As one can observe from figure 1.1 (i), the lymphatic network starts close to blood capillaries where fluid is pumped and releases lymph back to the venous system, after passing through few lymph nodes. In this system one can categorize into three main components:

• Lymph vessels (which are described below),

• Lymph,

• Lymph nodes.

Lymph vessels form a network that covers most of the human body. The skin, intestines and lungs are a few example of organs covered with the lymphatic network. Figure 1.1 (ii) represents the inside of a right leg skin collecting lymphatic network. As it can be seen, the bottom of the foot is covered with vessels which extend up to the lymph node of the groin area. There are also some hidden nodes under the knee called Popliteal lymph nodes. As defined by Sappey [1874], one can divide the network into three hierarchical levels corresponding to increasing diameter ranges:

• Capillaries (also called initial lymphatics) are smooth channels with no valves and an average diameter of 50 µm Sappey [1874].

• Secondary vessels (also called collecting lymphatics) have a diameter of 0.12 mm on average Sappey [1874] or 0.2 mm depending on authors [START_REF] Margaris | Modelling the lymphatic system: challenges and opportunities[END_REF]. They are made of units called lymphangions that can pump fluid thanks to the presence of valves, with an inter-valve distance/diameter ratio of about 10. These will be described in more detail below.

• Trunks which have the same structure as secondary vessels but with a lower inter-valve distance/diameter ratio. They often run parallel to the largest veins and are less numerous than the secondary vessels.

The focus of this thesis is on the secondary vessels, which are also called collecting lymphatics. This terminology has been used by [START_REF] Zawieja | Lymph protein concentration in initial and collecting lymphatics of the rat[END_REF] and more recently in [START_REF] Breslin | Lymphatic Vessel Network Structure and Physiology[END_REF]. Lymphangions are the elements constituting the lymphatic network and are represented in figures 1.3 and 1.1 (iii).

The figure also depicts the length of a lymphangion and its diameter. As they contract and expand to propel lymph which is contained inside, their diameter varies. The cycle used by the pump is described in the 0D model chapter. How-ever it is possible to get an average value around which the diameter will be considered contracted or relaxed. Table 1.1 shows data from the literature, in which the value of 0.2 mm is often used for lymphangion diameter. Furthermore the length depends on the source, but a value of 2 mm will be used and seems to be a consensus. It is important to note that human lymphangions are not easily accessible and requires experimental know-how since lymph vessels are smaller than blood ones and also present in smaller quantities. Furthermore, ethic measures are required to observe those vessels in vitro. However with technical progress, it is possible to observe them via specific markers but it is not still widely spread. In literature, lymphangions of rats or other mammals are often more studied than human ones. For example [START_REF] Arkill | The structure and mechanical properties of collecting lymphatic vessels: an investigation using multimodal nonlinear microscopy[END_REF] used bovine lymphatics, Bohlen et al. [2009] worked on rat lymphatic valve, or even dogs [START_REF] Burton-Opitz | The viscosity of lymph[END_REF]. The lymphatic network begins with a lot of porous dead-end capillaries which merge together to form a bigger lymphatic vessel. From then, the lymphatic vessels can merge and fork at lymphatic nodes, or without lymphatic nodes, as can be observed on figure 1.2. The walls of a lymphangion are composed of lymphatic smooth muscles according to von der [START_REF] Weid | Lymphatic smooth muscle: the motor unit of lymph drainage[END_REF]. The thickness in table 1.1 is defined as 10 µm, however as the wall contracts and relaxes, this value should be considered as a variable. These walls are essential for the drainage of interstitial fluid through the contraction cycle. Lymphangion walls behave as a feedback loop from a wide variety of stimuli that helps to modulate their contractile activity. As of today it is still unclear whether stimuli are received from chemical balance (Na/CO ions) or electro-nervous impulse.

Besides, from an energy point of view, only 40% of the fluid displacement is generated by these contraction cycles. The rest of the displacement is ensured by muscle contraction due to body movement or other internal displacements squeezing the lymphatic vessels. One can say that physical activity is essential for lymph circulation. Figure 1.3 shows a schematic of a lymphangion. The valve, which in our case is located on the left of a lymphangion is attached to the wall and floating in the lymph. It is formed of two leaflets, their shape like a moon croissant and means of attach to the lymphangion itself favors the flow in one specific direction, so one can say that these are unidirectional valves. According to [START_REF] Mazzoni | Structure of lymphatic valves in the spinotrapezius muscle of the rat[END_REF], there is no active smooth muscle action inside the leaflet tissues to open or close the valve, making it a purely passive object. Furthermore, only the pressure and the viscous forces appear to displace the valve, the valve is constructed such that it prevents prolapse and therefore reduces the risk of regurgitation (back-flow). Also in [START_REF] Daniel | Integrated geometric and mechanical analysis of an image-based lymphatic valve[END_REF], they were able to reconstruct a valve using a confocal microscope method for a rat lymphatic valve. It confirms that the valve is bicuspid and that each valve has a semilunar shape. This shape has a bias toward forward flow and tend to close in case of reverse flow. The valve is composed of elastic fibers and collagen matrix according to [START_REF] Arkill | The structure and mechanical properties of collecting lymphatic vessels: an investigation using multimodal nonlinear microscopy[END_REF]. In terms of elastic properties, the Young modulus is estimated to be around 45 kPa by [START_REF] Wilson | Determining the combined effect of the lymphatic valve leaflets and sinus on resistance to forward flow[END_REF], however [START_REF] Wilson | A fully coupled fluid-structure interaction model of the secondary lymphatic valve[END_REF] propose that it might be lower and should be around 20 kPa.

Now that the network and secondary lymphatic structure have been defined, the next section is about lymph, the fluid itself. So it will expose its composition properties, as well as the function it covers as a fluid.

Lymph and its composition

Lymph is a fluid (lympha meaning water, and also related to water in movement in Latin) issued from the filtration of capillaries which are the dead-end entrance point of the network. As mentioned by [START_REF] Moore | Lymphatic system flows[END_REF] the total flow rate for a human person should be around 5 L day -1 , and the total amount of lymph is between 8 to 10 liters according to [START_REF] Moore | Clinically oriented anatomy[END_REF]. As a comparison, blood flow rates vary between 3.0∼26mL min -1 in arteries and 1.2∼4.8mL min -1 in veins according to [START_REF] Klarhöfer | High-resolution blood flow velocity measurements in the human finger[END_REF] Furthermore the total blood flow rate is 5 L min -1 . The composition of lymph evolves while crossing lymph nodes. As a matter of fact, pre-nodal lymph is similar to interstitial fluid and its composition is mainly water, salts, plasma proteins and white blood cells. The concentration for plasma proteins and cells is so low that one can consider lymph to be as viscous as water in that case [START_REF] Stoltz | Letter: Study of the rheological properties of human lymph[END_REF]. The viscosity of lymph has been found to be 1.7 times greater than water according to [START_REF] Burton-Opitz | The viscosity of lymph[END_REF]. Before flowing into the capillaries, the fluid is called interstitial fluid, which allows interstitial cells to interact with hormones or proteins Fogh-Andersen et al. [1995]. The osmotic pressure difference plays a key role in the feeding of lymphatic capillaries by the interstitial fluid Guyton et al. [1976]. However, regarding the size of opening of the capillaries, some molecules such as glycocalyx play a role of sieve in the formation of lymph according to [START_REF] Kirk | Lymph formation, composition and circulation: a proteomics perspective[END_REF]. As long as interstitial fluid hydro-static pressure is above lymph pressure, this pressure difference allows interstitial fluid to be drawn inside the lymph capillaries. Furthermore, lymph can be seen as issued from the filtration of blood to the interstitial fluid as mentionned in [START_REF] Quéré | Description anatomique et histologique, physiologie du système lymphatique[END_REF]. The fluid enters the lymphatic capillaries with a flow rate of approximately 5 mL h -1 as mention in [START_REF] Kirk | Lymph formation, composition and circulation: a proteomics perspective[END_REF] and [START_REF] Smith | The traffic of cells through tissues: a study of peripheral lymph in sheep[END_REF], from a vessel entering a lymph node. However, this value is dynamic and it has been observed to vary, notably under pathological conditions. Those variations can be associated with increased lymphangiogenesis, cellular trafficking or inflammatory reaction, which affect the lymphatic system.

From a functional viewpoint, lymph is considered as part of the immune system. As mentioned in [START_REF] Avi | Lymphocyte trafficking[END_REF], lymph is charged with lymphocytes and monocytes. Between the capillaries and the collecting system, the concentration of these cells in lymph can increase if the immune system reaches inflammatory conditions. In addition, the stiffness of lymphocytes can evolve as it is discriminated toward a specific purpose, as shown in figure 1. 4. After going through a lymph node, more antigen/antibody complexes can be found in the lymph.

Also surrounding the gut and during the digestion process, the lymphatic system carries chyle. According to [START_REF] Joseph | On the composition of chyle chylomicrons[END_REF], this fluid is specifically produced at the wall of small intestine and differs from lymph because it contains emulsified fats. These lipids are stored in chylomicrons, which form particles in the fluid forming a suspension. As the concentration of fat is greater than in normal lymph, chyle is whiter and opaque. Furthermore the viscosity of chyle is expected to be different from lymph and some study on drug absorption hinted the viscoelastic behaviour of such fluid [START_REF] Peppas | A theory of molecular diffusion in the intestinal mucus[END_REF]. Another point to raise in the modification of lymph is the alteration of lymphocytes. As they are part of the immune system, for this section only the B and T cells are considered. From a general point of view, B lymphocytes play a large role in the adaptive immune response, while T-cells are more active during the primary response, and can also activate B-cells, as explained in DeFranco [1987]. Furthermore, it has been shown by [START_REF] Learmonth Gowans | The route of re-circulation of lymphocytes in the rat[END_REF] that lymphocytes primarily reside in lymph nodes. They massively enter the blood circulation and lymph circulation if the immune system is provoked by the intrusion of pathogen agents. There is information on density of lymphocite in rat lymph, the density ranges from 32.6 lymphocites per µL to 3.55 × 10 4 according to [START_REF] Dixon | Lymph flow, shear stress, and lymphocyte velocity in rat mesenteric prenodal lymphatics[END_REF]. No information on human lymphocites density in lymph has been found and it is expected to vary between rat and human or also depends on the health between different specimens. Lymph nodes also allow for precise and localized reaction by connecting lymph and the blood circulation.

As shown in figure 1.1, the lymphatic system collects fluid, cleans it and inserts it back to the main blood circulation. This part of the network can fail and lead to multiple disease, such as lymphoedema. Complications from this disease can become significant if patients are sedentary Greene [2015]. In terms of treatment, massage, inelastic compression bandage and skin care are a few examples of interventions in order to manage the complications Ko et al. [1998]. Other surgeons have investigated rerouting the network in order to achieve drainage of lymph after removal of a lymph node (lymphadenectomy). This operation often happens in breast cancer, and can leave some lymphatic vessels open which might then lead to lymphoedema. However, there are no clear standard or efficient surgical techniques for such operation to this day [START_REF] Scaglioni | Systematic review of lymphovenous anastomosis (lva) for the treatment of lymphedema[END_REF].

Outline of the Thesis

The main goal of this thesis is to gather knowledge on the physical phenomena that occur within the lymph network. As of today, it is not possible to cure a lymphoedema, and physician often use compression garments to prevent this long-term condition. So by improving the knowledge of lymph transport, one should find better ways to prevent or even cure lymphoedeme.

From a physical point of view, the fluid and the solid particles inside lymph or the flexible walls of the vessel present new difficulties and unsolved problems when trying to model and understand the transport of lymph. Some key challenges for which this work seeks to find some further insights are:

• The pumping function of lymphangions and its dependency on fluid properties, mechanical properties of the lymphangion and valves and dynamical parameters such as pressures and contractility

• The structure of the flow of a particle-laden fluid (lymph and lymphocytes) in lymphangions: how does the presence of particles (cells) affect the structure of the flow, its rheological properties and the pumping efficiency

• The modeling of lymphangion networks viewed as dynamical systems of coupled active units, their ability to synchronize to produce a global flow and the description of flows in such complex networks

In order to contribute to these different aspects, two different scales of the system were considered with two types of modeling, namely 2D fluid-structure simulations at the scale of a lymphangion, and a 0D model of networks in which each lymphangion is modeled by a set of coupled ordinary differential equations.

To conclude this introduction, the main chapters of this thesis are listed hereafter. The first chapter provides a review and gathers multiple articles with model or modeling made similar to the work presented in this thesis. Then, the second chapter will develop an improved version of the model in 0D dimension derived by [START_REF] Bertram | Simulation of a chain of collapsible contracting lymphangions with progressive valve closure[END_REF] and its use against an example of a leg skin network. The third chapter focuses on the description of a fluidstructure interaction model and 2D simulation code that was used in the fourth chapter. The last chapter presents a numerical study on the influence of the presence of lymphocytes on lymph flow in a vessel, with a focus on the dispersion and transport of particles through a lymphangion with pulsed flow. Finally, conclusions and perspectives for future work are presented.

Introduction

The human lymphatic system can be described of three main components, which are the lymph, the network and the lymph nodes. Lymph is a body fluid, mostly made of blood plasma, its molecular composition is subject to many variations as it helps to spread hormones D' Alessio et al. [2007] or also to support the immune system [START_REF] Olszewski | Lymphatics, lymph and lymphoid cells: An integrated immune system[END_REF]. Furthermore, one of the most notable case is the modification in lipid concentration within the lymph near the guts after the digestive process. The lymphatic network can be separated into three main categories: the capillaries, the vessels and the trunks. The former is valveless and quite porous. The latter two have the same structure, a succession of lymphangions and only differ by their diameters. They both are able to contract due to muscle cells present in their walls. Lymph nodes are of oval shape, they 13 play key role in the immune system as they produce (as well as bone marrow) and conserve different version of antibodies Jerne [1973].

For this thesis, we will focus on the mechanical study of the lymphangions, its pumping mechanism at the lymphangion scale and into the network. Lymph gets into the network from the interstitial fluid into the capillaries. In fact, using osmotic pressure the interstitial fluid flows into the porous end of the lymphatic capillaries. According to [START_REF] Philibert | Anatomie, physiologie, pathologie des vaisseaux lymphatiques considérés chez l'homme et les vertébrés[END_REF] the diameter of those are lower than blood capillaries. For most of them the diameter is not larger than a micrometer. All these capillaries converge to larger vessels which are regularly valved and have muscle tissue. These valves are important for the flow direction. In fact, they close to prevent back flow and open when the flow pushes them open, in other words, these valves have a passive bias for forward flow. Furthermore, there are two types of forces to generate lymph transport. Passive ones are generated by overall movement of the entire body such are walking, exercising. It can be linked to the contraction of skeletal muscles, fluctuation of the central venous system, the influence of the respiratory movements, gravitational forces Gashev [2002] or more simply by walking or exercising. Active forces are generated by the muscles in the lymphatic vessels or trunks. Self-contractile activity has been observed in the lymphatic vascular system.

In order to get overall value of the forces per lymphangion or global scale of the total flow of lymph without measuring it in vitro, which might disturb the measurement, a zero dimension model can be easily set up in order to simulate a whole collecting network. Different models are shown in detail then different papers using networks is presented. Finally, the lymphangion approach with two or three dimensions using fluid structure interaction simulation are shown.

Various lumped models

The initial work on lymphatic lumped modelling is from [START_REF] Narender | A computer model of the lymphatic system[END_REF], and used main lymphatic trunks to form a network. However, trunks and secondary vessels are different in length-diameter ratio. The main equations are common to most of the articles since each model the same phenomenon. These equations can be written as:

dV dt = Q in -Q out Q = -Kδp p = p act + p pas + p ext (2.1)
This set of equations consists in mass conservation, the Hagen-Poiseuille law and the pressure balance involving fluid pressure p active pressure p act (wall contraction), passive pressure p pas (wall elasticity) and external pressure p ext . With p the pressure in the lymphangion, V the lymphangion's volume, Q the flow and K i the conductance (resistance's inverse) with the subscript i as the valve or the lymphatic, since those have independent diameters. A complete derivation of this set of equations is provided in the following chapter.

In [START_REF] Bertram | Simulation of a chain of collapsible contracting lymphangions with progressive valve closure[END_REF] from a modelling point of view, a simple geometry is defined with the wall of a lymphangion forming a cylinder. The pressure drop is estimated from the standard Poiseuille flow in a straight channel. Such approximation is strong since the lymphatic is displaced via a pulsed flow. However the time spent during the contraction (systole) and the relaxation (diastole) are not equal. For lymphangions, the diastole is significantly longer than the systole in a normal condition [START_REF] Meisner | Lymphatic vessels transition to state of summation above a critical contraction frequency[END_REF], [START_REF] Telinius | Human thoracic duct in vitro: diameter-tension properties, spontaneous and evoked contractile activity[END_REF]. So the lymph is often in a steady state, thus the Poiseuille flow approximation is valid. The valve is defined in a similar manner, the pressure difference between the two ends of the valve determines the valve resistance, then a pressure drop is computed. To sum up, this method considers a lymphangion as a chain of two cylinders of different diameter with flow passing through. Then, in order to compute pressures this paper defines the transmural pressure. To simulate the contraction and relaxation process of the lymphangion's wall, the constitutive equation has a forced cycle with a trigonometric function. Furthermore, the relaxation term is the sum of an exponential term and polynomial series to fit the pressure-diameter relation given by the experimental data, and the collapse function is using an inverse cubic term. To solve the set of equations, a Runge-Kutta fourth-order method is used for time integration and Matlab own subroutine for ordinary differential equation (ODE) is used. [START_REF] Contarino | A one-dimensional mathematical model of collecting lymphatics coupled with an electro-fluid-mechanical contraction model and valve dynamics[END_REF], use a one dimension approach, a set of partial differential equations (PDE) similar to (2.1) however its discretization does not depend on the lymphangion. In other words, this approach does not hold the variable value at specific points of the lymphangions, such as before or after the valve. Also, this paper introduces a model for the electro-fluid-mechanical contraction, which takes into account amount Ca 2+ in lymph to increase rate of contraction while NO help modulate such response and play the role of contraction inhibitor1 . This allows the contraction term in the constitutive equation to be independent of any trigonometric periodic function referred by s(x, t) (the contraction state variable) in table 2.1. Instead the cycle is based on an equilibrium between the stretch-activated calcium influx and contraction inhibitor of the wall shear stress by the production of nitric oxide. This equilibrium provides its own periodic behaviour, so the lymphangion contracts and relax depending on its input and output pressure. It is solved using an ordinary differential first order in time equation with s(x, t) (a contraction state) as unknown, explained in table 2.1. Tests are run on a chain of three lymphangions in a similar set up to [START_REF] Bertram | Simulation of a chain of collapsible contracting lymphangions with progressive valve closure[END_REF] paper. Furthermore, the valve model is implemented in a similar manner. A sensitivity analysis is done for different pressure gradients, and the presence of electrical coupling in order to pass contraction information to synchronise lymphatic contractile activity is suggested in the conclusion. [START_REF] Jamalian | Parameter sensitivity analysis of a lumped-parameter model of a chain of lymphangions in series[END_REF], uses a model similar to [START_REF] Bertram | Simulation of a chain of collapsible contracting lymphangions with progressive valve closure[END_REF] which assumes Poiseuille flow in lymphangions. Furthermore, the contraction term of the constitutive equation is forced by a trigonometric function and does not include the balance between the wall shear stress and the contraction inhibitor to compute the contraction-relaxation cycles. Matlab subroutines are used to solve the system of equations. A sensitivity analysis is performed for the valve parameters and for lymphangion wall as well. In order to estimate the most efficient chain of lymphangions, their number in a chain is varied. Also, [START_REF] Jamalian | Parameter sensitivity analysis of a lumped-parameter model of a chain of lymphangions in series[END_REF] points out that rat lymphatic secondary vessels of human and rat should behave differently. The constitutive equation used in this article is written in table 2.1. The eleven parameters on the passive term, and the four parameters for the active term makes it the most complex constitutive equation of the literature. [START_REF] Tretyakova | Modelling lymph flow in the lymphatic system: from 0d to 1d spatial resolution[END_REF] defines two similar models, the first in zero dimension and the second in one dimension. For the zero dimension model, the constitutive equation used in this article is mentioned in table 2.1. Both of them use a similar constitutive contraction equation such as [START_REF] Bertram | Simulation of a chain of collapsible contracting lymphangions with progressive valve closure[END_REF]. The valve modelling however varies quite a lot from the other papers. In fact, the model is implemented in order to prevent back flow on a general level, and does not have a dynamical valve model. The zero dimension model has a constant value as a Poiseuille resistance, the one dimension model applies friction force to compensate back flow. The geometry used for the zero dimension is a series of lymphangions, from two to ten in a straight channel. In the description of the one dimension method, a network approach is mentioned, but the numerical section does present any results. The metrics used in this paper are pressure, radius flow and contraction term. [START_REF] Bertram | Simulation of a chain of collapsible contracting lymphangions with progressive valve closure[END_REF], [START_REF] Contarino | A one-dimensional mathematical model of collecting lymphatics coupled with an electro-fluid-mechanical contraction model and valve dynamics[END_REF], [START_REF] Jamalian | Parameter sensitivity analysis of a lumped-parameter model of a chain of lymphangions in series[END_REF][START_REF] Tretyakova | Modelling lymph flow in the lymphatic system: from 0d to 1d spatial resolution[END_REF] lymphangion models are often applied to small or larger networks in order to understand the dynamic of lymph flow. Simulations can help clinicians as they provide insight of flow and pressure of the lymphatic system without the need of experimental protocol or surgery operations.

Lymphatic network lumped models

There are few articles that focus on simulating various kinds of lymphatic network. As detailed in the introduction there are different types of network. Unlike the articles from the previous section, the geometry presented here has convergence and divergence of lymphatic vessels. For example, it is the case in [START_REF] Jamalian | Network scale modeling of lymph transport and its effective pumping parameters[END_REF] and [START_REF] Narender | A computer model of the lymphatic system[END_REF] papers. The first one defines a geometry made of chains of three vessels merging or diverging from the secondary lymphatic network. The second one has a geometry based on human lymphatic trunks. Such modification in the structure has various impact on the modelling and the simulation. But those can be summed up as how to handle the bifurcation and confluence of vessels. From a model perspective, these connections are handled by equating the pressure and by satisfying the mass conservation at the connecting lymphangion. The boundary conditions on [START_REF] Jamalian | Network scale modeling of lymph transport and its effective pumping parameters[END_REF] paper are set so the inlet pressure remains the same and is defined at 6 cm H 2 O. Furthermore, the outlet pressure is varied in different simulations and extends from -3 to 41 cm H 2 O. However, networks used in this article are

P = K(x, t)Φ(x, t) + Pext K(x, t) = 1 12(1 -ν 2 ) (Emin + s(x, t)(Emax -Emin))( h0 r0 ) 3 Φ(x, t) = ( A(x, t) A0(x) ) m -( A(x, t) A0(x, t) ) n + c( A(x, t) A0(x) -1)
A is the crosssectional area, only made of convergent vessels with a total of 3 connections, with four inlets and one outlets, as shown in figure 2.1. [START_REF] Jamalian | Network scale modeling of lymph transport and its effective pumping parameters[END_REF] paper, here it is drawn with three lymphangions per vessels. In the article, the amount of lymphangions per vessel is varied. Furthermore, the overall structure remains the same, such that the number of connection between vessels remains at three with four inlets and one outlet.

P = fp(D) + fa(D, t) + Pext fp(D) = P d [c1( Di c9 -c2) 2 + c3 exp(c4( Di c9 -c5)) + c6 + c7( Di c9 -c8) + c10( Di c9 ) -3 ] fa(D, t) = M0 Di (1 -cos(2πf (t -tc)))( 1 1 + exp(-s d (Di -Da)) + 1 1 + exp(s d (Di -Da)) -1)
C 1 C 2 C 3 Figure 2.1: Schematic of a network used in
In this article, [START_REF] Narender | A computer model of the lymphatic system[END_REF] which is the oldest one related to simulating lymphatic network, it does consider, the whole body main ducts going back to blood, and does not contain bifurcation. The network modeled is slightly less than 300 lymphangions. Pressure profile against time is observed for different organ. In order to simplify the network, all lymphangions are considered the same, which is used as well for our model. With regards to the figure provided in this article, the time axis shows ten seconds of simulation time. It is acceptable to consider this time not enough to propagate information into the entire network. Also the number of lymphangions can be considered small, since trunks ones have small ratio than secondary vessel ones.

Numerical method, fluid-structure interaction

In this section, we brielfy review the main numerical methods that can be used for numerical simulation of lymph flow. In order to solve numerically the displacement of lymph in two or three dimensions, one needs to solve the fluid and the solid displacement. There are multiple ways to couple fluid and solid solver, the mains approaches are:

• Arbitrary Lagrangian Eulerian (ALE),

• Smoothed Hydrodynamic Particles (SPH),

• Immersed Structural Potential Method (ISPM).

The ALE method finds its root with the work of [START_REF] Hirt | An arbitrary lagrangian-eulerian computing method for all flow speeds[END_REF]. It is limited since fluid and solid are meshed on the same space and need to be re-meshed when the solid is displaced, and can be computationally expensive [START_REF] Gil | The immersed structural potential method for haemodynamic applications[END_REF]. Meanwhile the SPH method does not have a mesh, each particles need to compute its neighbor distance in order to estimate the next time-step which is demanding. However, this process is highly efficient with proper parallelization. The immersed boundary method which founded the ISPM was first introduced by [START_REF] Charles S Peskin | Flow patterns around heart valves: A numerical method[END_REF]. ISPM is composed of a mesh for the fluid and a cloud of particles for the solid. This means that for every time-step the information of velocity needs to be interpolated from the fluid grid to the solid points raising the computational time. A comparison on FSI benchmark, a dam break with elastic obstacle has been done by [START_REF] Yang | Unified one-fluid formulation for incompressible flexible solids and multiphase flows: Application to hydrodynamics using the immersed structural potential method (ISPM)[END_REF], of between SPH and ISPM.

These methods found an echo in various domains from aero-elasticity to mechanical biology. However, the time scales varies a lot between these two domains. Naval applications is confronted to similar problem of elasticity as aerodynamicist. In fact, the Reynolds number for ships or planes is larger by three up to five order of magnitude compared to biological problems, thus inertia plays a larger role. As detailled in [START_REF] Sigrist | Fluid-Structure Interaction[END_REF], fsi coupling can be written in the Laplace space. This frequency based formulation is suited for aerodynamic elasticity problems, as the solid does not deform a lot but with higher frequency than in biological application.

In the mechanical biology domain, few problem has been studied such as the aortic valve or red blood cells. The first replacement of this valve was first performed in 1962. Since then, teams try to improve the design, but the in vivo measurement and working conditions are difficult to obtain valuable results. [START_REF] Thubrikar | The design of the normal aortic valve[END_REF] presented work on improving the working functionality of an aortic valve on dogs. But with the help of simulation and the increase of computational power available, not only valve were considered but also blood vessels. The figure 2.2 (i) & (ii) represent an aortic valve from a numerical point of view Such that today, there are frameworks to compute intraventricular blood flow [START_REF] Prud'homme | Feel++ : A computational framework for galerkin methods and advanced numerical methods[END_REF] or [START_REF] Spühler | 3d fluid-structure interaction simulation of aortic valves using a unified continuum ALE FEM model[END_REF]. However, on the lymphatic side progress are moderate, the secondary network valve has been observed using confocal imagery and mechanically analysed by [START_REF] Daniel | Integrated geometric and mechanical analysis of an image-based lymphatic valve[END_REF]. Also, on the next section we are investigating a new approach to the lymphatic valve problem.

Two & three dimensions, secondary lymphatic valve and lymphangion models

In this section, we are going through articles focusing on the numerical simulations of lymphangion using coupled or uncoupled fluid-solid interaction (FSI). By opposition to the lumped model simulations, two and three dimensions simulations are more recent in the literature. Most of articles on this topic focus on the valve displacement and its behaviour in the pumping function of a lymphangion. They also use water as a fluid in the simulations, by setting the density and the viscosity, apart from [START_REF] Bertram | Modelling secondary lymphatic valves with a flexible vessel wall: how geometry and material properties combine to provide function[END_REF]. [START_REF] Wilson | Determining the combined effect of the lymphatic valve leaflets and sinus on resistance to forward flow[END_REF] uses an uncoupled approach as forces are applied to the valve structure giving a deformation. This deformation uses then as input geometry for a fluid simulation. The limitation of such method to compute the valve resistance has clear limits by not coupling fluid and solid. However, the geometry are reconstituted from confocal microscopy images. The valve stiffness is homogeneous in the valve and is set to 45 kPa. The stiffness of the valve used here is relatively high compared to other paper as shown in table 2.3.

(i) (iii) (ii)
It focuses on the pressure along side the central axis of the lymphangion, which exhibits a typical Poiseuille flow in a straight tube type of pressure profile with variation near the valve. A later paper from the same group, [START_REF] Wilson | A fully coupled fluid-structure interaction model of the secondary lymphatic valve[END_REF] achieved a fully coupled transient simulation of flow through a valve. While remaining focused on the sinus region, it presents the fluid dynamic, valve deflection and the resistance to forward flow. However, the geometry used in the simulation only corresponds to one forth of the channel, based on its symmetrical nature, which might not completely capture the physic, since the Reynolds number is low and eddy asymmetry will not establish. Also, vessel walls are not simulated, yet they are likely to influence the dynamic of the valve.

Using a different method, but with a very similar geometry, [START_REF] Ballard | Probing the effect of morphology on lymphatic valve dynamic function[END_REF] focuses on impact of valve aspect ratio and the relative stiffness of the valve. In fact, it uses a lattice Boltzmann method for the fluid is coupled to a lattice spring model for the solid. This method is described in [START_REF] Ladd | Lattice-boltzmann simulations of particle-fluid suspensions[END_REF]. From a geometry perspective, the walls are statics and only the valve is moving. Furthermore, [START_REF] Ballard | Probing the effect of morphology on lymphatic valve dynamic function[END_REF] mention observing valve with aspect ratio from 1.15 to 3.0, and [START_REF] Wilson | Determining the combined effect of the lymphatic valve leaflets and sinus on resistance to forward flow[END_REF] have aspect ratio from 1.17 to 1.41. Their valve stiffness is homogeneous but varies from simulations from 1.2 kPa to 7.5 kPa. By opposition to the first paper, [START_REF] Ballard | Probing the effect of morphology on lymphatic valve dynamic function[END_REF] use a different valve Poisson ratio of one third, while all the other papers use one half in term of ratio. Aside from the simulation, experimental work on lymphatic vessels from rats and sheep are presented. [START_REF] Bertram | Modelling secondary lymphatic valves with a flexible vessel wall: how geometry and material properties combine to provide function[END_REF] and [START_REF] Li | The effects of valve leaflet mechanics on lymphatic pumping assessed using numerical simulations[END_REF] work include computation of wall deformation. This brings new insight as the vessel wall is responsible for part of the lymph drainage since it produces force to squeeze the fluid out. [START_REF] Li | The effects of valve leaflet mechanics on lymphatic pumping assessed using numerical simulations[END_REF] use lattice Boltzmann model for the fluid, furthermore they track ions couple Ca 2+ and NO to trigger respectively relaxations and contractions of the vessel walls. The NO concentration is treated as a scalar field, it uses convective, diffusive and an additional term to take into account the decay and production rate, shown in the equation number 7. Furthermore, the Ca 2+ is only confined to the wall of muscle cells, and only transfer from one cell to its neighbors, thus the concentration equation is only diffusive. It is more complex than the nitrite oxide equation, due to its non linearity and fail-safe feature to limit the minimum of radius. From the two dimensions geometry, valves need variation of Young's modulus to correct the parabolic anchoring to the wall which can only be properly defined using three dimensions geometry. Here, the valve stiffness at the trailing edge is ten-times smaller than at the root, and sinus region in the result's figures does not exhibit the round shape typical of lymphatic vessels. It is not clear whether these two phenomena are linked. [START_REF] Li | The effects of valve leaflet mechanics on lymphatic pumping assessed using numerical simulations[END_REF] shows improvement of pumping efficiency with valve trailing edge flexibility, and that the valve rest position does impact the pumping. [START_REF] Bertram | Modelling secondary lymphatic valves with a flexible vessel wall: how geometry and material properties combine to provide function[END_REF] runs its simulation with Arbitrary Lagrangian-Eulerian approach, using a three dimensions geometry. However, it is not clear if the authors used a full model of a lymphangion or 1 4 of the whole geometry. It shows evidence of retrograde flow during the valve closure, and rarely closes completely. Furthermore, the vessel wall passively deforms and does not account for the muscular forces, which can impact the pressure inside the lymphangion and thus the valve shape. Since the valve is incompressible and its deformation due to pressure difference for it to open or close, it must have an impact on the wall displacement as well. Valve resistance, flow rate and valve deformation are computed against pressure difference in order to compare with lumped models and verify their accuracy.

Method

Geometry Metrics [START_REF] Bertram | Modelling secondary lymphatic valves with a flexible vessel wall: how geometry and material properties combine to provide function[END_REF] Adina software using Arbitrary Lagrangian-Eulerian (ALE) approach 3D Both wall and valve are simulated. Variation in the wall and valve stiffness are achieved.

• Valve trailing edge Abstract -We present a flow model of the lymphatic network based on a zero dimension model of the lymphangion dynamics with a new numerical formulation, for calculation of lymph flow in the collecting network. This study starts on a simplified version of the model that reduces the number of parameters in the constitutive equation without altering its physics. It will focus on different specific cases such as confluences or bifurcation geometries which are key elements necessary for the creation of a network. The results will be compared with experimental data acheived on laboratory rats. Finally, the results of a specific network geometry are presented from an anatomical drawing of a human leg.

Introduction

In this chapter, a 0D model is presented for the estimation of lymph flow in chains of lymphangions which is then used to study flows in small-scale networks as well as in a larger network whose geometry was extracted from anatomical drawing from Sappey [1874] This model is an improvement of the previous one presented in [START_REF] Jamalian | Network scale modeling of lymph transport and its effective pumping parameters[END_REF] as it does not enforce the pulsatile behaviour 25 by trigonometric function. We will use it against experimental data as a validation procedure. Most of the experimental data used in this thesis are produced either in [START_REF] Michael | Intrinsic increase in lymphangion muscle contractility in response to elevated afterload[END_REF] or [START_REF] Bertram | Consequences of intravascular lymphatic valve properties: a study 135 of contraction timing in a multi-lymphangion model[END_REF]. It focuses on the main lymphatic channels, which Sappey [1874] calls trunks. They differ from collecting lymphatics by their length and diameter. In fact, when new lymphangion mechanisms are observed by experimentalists, new models are developed to try to fit those experimental data. As well as a vessel made of a chain of a few lymphangions (most of the time ten) was studied experimentally by [START_REF] Bertram | Incorporating measured valve properties into a numerical model of a lymphatic vessel[END_REF] and numerically by [START_REF] Bertram | An improved model of an actively contracting lymphatic vessel composed of several lymphangions: pumping characteristics[END_REF] or [START_REF] Tretyakova | Modelling lymph flow in the lymphatic system: from 0d to 1d spatial resolution[END_REF]. Note that, gathering mechanical data from real lymphangions requires fine equipment with proper tuning in terms of size, forces and chemical balance. The model used in this section allows computer simulations of a lymphatic network using a 0D approach. It is based on a few approximations such as:

• lymph being a Newtonian fluid with the viscosity of water,

• flow being in a quasi-steady state, in other words, the flow is assumed to be slowly varying, such that a fully developed Stokes flow is assumed,

• lymphangion geometry is considered to be cylindrical.

In terms of dimensionless numbers the Reynolds, the Womersley and the Strouhal numbers are defined as:

Re = vDρ µ , Wo = D ωρ µ and, St = ωD v (3.1)
With characteristic length D = 0.3 cm the diameter of the vessel, the density ρ = 1 g cm -3 the dynamic viscosity µ = 1 cP and ω = 0.01 Hz for pulsation frequency. All these values are gathered from the more exhaustive table 3.3. The frequency used here fit a human lymphatic system but in this document a comparaison is done with experimental data which uses rats or mouses.

As the proposed code does not compute the velocity but pressure, it is possible to use the pressure gradient over a Poiseuille flow in a straight channel the characteristic velocity v = 1.36 cm s -1 , using the following equation:

v = D 2 8µ δp δx
With the pressure gradient being the difference of pressure boundary condition over three lymphangions length. Using these values one can determine the Reynolds number to be 40 and the Wormersley number 0.3 approximately. The Reynolds number correspond to the ratio between the inertial forces to the viscous forces, while the Womersley number compares the ratio of inertial pulsatile flow effects to viscous effects. So here, the Reynolds number is low as well as the Womersley number. Therefore, the unsteadiness of the pulsatile flow is not considered as the Reynolds number is low. Moreover, when Wo is small, the velocity profile solution is similar the Poiseuille profile. The Strouhal number compares the ratio of velocity between oscillatory velocity to ambient flow velocity, it can be written as a combination of both the Reynolds and the Womersley number. For convenience, and because most of the literature use CGS (centimeter-gram-second) unit system, this unit system will be used in this document. Also, it is important to notice that the pressure and work will be considered in this CGS system, if not indicated otherwise. The table 3.1 introduces the physical quantities, their CGS units and their SI equivalents.

The next section is about the presentation of the model. Then we will see how this model is implemented. Finally, the results section will discuss different cases run using this model.

Description of a lymphangion

A lymphatic vessel is made out of a chain of lymphangions. A single lymphangion can be simplified as shown in figure 3.1. In all points described in figure 3.1, the [i] subscript represents the i th lymphangion.

i th lymphangion i -1 th i + 1 th P n[i-1] P v[i] P [i] P n[i] D [i] D v[i+1] Q 1 Q 2 Q 3 (i)
Volume V wall motion v[i+1] are represented only to show the iteration over the next computational unit that leads to lymphangions if put together sequentially. In order to properly describe how the lymphangion is defined, first the modelling equations for the valve are presented, then the equations governing the contraction of the lymphangion. From these two models, the coupling between these two sets of equations is explained.

Valve model

The model used here, is based on the articles by [START_REF] Mynard | A simple, versatile valve model for use in lumped parameter and one-dimensional cardiovascular models[END_REF] & [START_REF] Bertram | Simulation of a chain of collapsible contracting lymphangions with progressive valve closure[END_REF] and is utilized to determine the diameter of the valve orifice. This computation is based on the pressure difference across the valve and the diameter of both the lymphangion and the valve interface. The valve state ζ is defined as: ζ ∈ [0 , 1] where 0 represents a fully open valve and 1 represents a fully closed one. The following convention is set to be 0 for a fully open valve and 1 for a fully closed one. ζ depends on time and the pressure difference across the valve. In order to compute the valve state, one needs to know whether it is closing or opening. This information is obtained by computing the difference of pressure across the valve. The pressure computational points is made of three different

P n[i-1] P v[i]
k open opening compliance (Ba -1 s -1 ) k close closing compliance (Ba -1 s -1 ) Figure 3.2: Sketch of a lymphangion's valve, with the mesh pressure point. points for each lymphangion: P v[i] , P [i] and P n[i] (cf. Figure 3.1 and Figure 3.2). Thus, the difference in pressure is computed as :

δP = P v[i] -P n[i-1] ,
with the subscript indicating the previous ( [i-1] ) or the current ( [i] ) connected lymphangion. Hence, if δP is positive the valve should be closing and negative the valve should be opening. Furthermore, the valve state is time dependent, one can define its derivative: dζ dt . If the pressure difference is positive, one can observe the valve is closing, and if it is negative, the valve is opening. Then, using the simplified model defined in the article from [START_REF] Mynard | A simple, versatile valve model for use in lumped parameter and one-dimensional cardiovascular models[END_REF], one has the following set of equation. Then in order to obtain the valve orifice diameter (D v ), the following equation is used. Here, D [i-1] and D [i] refers to the diameters of the two lymphangions surrounding the valve.

dζ dt =      (1 -ζ)k close δP if (δP ) > 0,
D v[i] = 1 2 (D [i-1] + D [i] )(1 -ζ) (3.3)
The previous equation requires the orifice diameter to be lower or at most equal to the average size of its two closest lymphangions. In most of the simulations, k close and k open hold the same value and is set to 100 Ba -1 s -1 . The characteristic time for opening and closing are therefore relatively similar. However, these also depend on the amplitude of δP , which depends on different factors, such as the state of contraction of preceding and following lymphangions. It is interesting to compare the contraction/relaxation cycle period with the valve characteristic time τ 0 .

τ 0 (t) = 1 kopenδP (t)
τ 0 characterizes the response of the valve to a given constant pressure difference input. In terms of scale, a human lymphangion contracts on average once every 20 seconds for an intestinal lymphatic vessel according to [START_REF] Telinius | Spontaneous and evoked contractility of human intestinal lymphatic vessels[END_REF].

Contraction model

Following the same logic as the valve state, equation (3.4) is the contraction intensity equation. ξ holds the value of the state of the contraction, is bounded between e min m 0 and e max m 0 , and is used in the active pressure formulation P a in equation (3.10). Detailed definition will be given in the following paragraph. γ is a variable used to contain ξ between two extrema and its definition is detailed later in this section.

ith lymphangion

P [i] D [i]
ω frequency (Hz) ξ contraction intensity (no unit) γ contraction memory (no unit) The contraction model used here aims to replicate the experimentally observed mechanism, which correlates the frequency of contraction ω to the diameter D [i] . There is a relation between Ca 2+ and NO and the contraction compliance of a lymphangion which explains this correlation, as stated by [START_REF] Contarino | A one-dimensional mathematical model of collecting lymphatics coupled with an electro-fluid-mechanical contraction model and valve dynamics[END_REF]. Such that Ca 2+ increases the rate of contraction while NO modulates it as a contraction inhibitor. However, in most lymphatic lumped models, the contraction is constrained in time by a periodic function (cf. [START_REF] Jamalian | Network scale modeling of lymph transport and its effective pumping parameters[END_REF], [START_REF] Tretyakova | Modelling lymph flow in the lymphatic system: from 0d to 1d spatial resolution[END_REF]). Here, the contraction is governed by an equation similar to the valve state equation.

dξ dt = (1 -ξ) k c ω if γ = 0, -ξ k r ω if γ = 1. (3.4)
ω is the frequency of contraction, it is linked to the diameter of the lymphangion as shown in equation 3.5, with c 1 the threshold for D D d to increase frequency and c 2 the frequency increasing rate.

ω = ω min + D D d -c 1 * c 2 if D D d ≥ c 1 (3.5)
The difference between this model and the valve model comes from the condition on γ. In fact, this condition relies on the comparison of diameter constant (cf. table 3.2) to the previous contraction intensity. γ is defined by equation (3.8). It can either have the value of 0 or 1, like a boolean and the process in determining its value will be described in the following paragraph. First, the contraction compliance k c and the relaxation compliance k r are defined as: e min is the minimum stiffness and e max is the maximum stiffness. e cen is the eccentricity of contraction and this parameter is bounded between 0 and 1. It defines the ratio between the time spent contracting and relaxing. m 0 is a stiffness factor for the active pressure. Standard values of those constants are given in table 3.3 From the differential equation (3.4), γ = 0 defines a contracting state while γ = 1 refers to a relaxing state. Such that, the contraction intensity ξ remains contracting if it has not reached e max m 0 , and it will remain relaxing while e min m 0 is not reached. It can be defined as follows:

k c = - e cen +
γ [i] =            0 if e min m 0 > ξ, 1 if e max m 0 < ξ, γ [i-1] else. (3.8)
This makes the lymphangion behave in a contractile fashion which changes whenever γ is above or below a trigger value. Equation (3.8) ensures that while the contraction memory(γ) remains constant, the contraction intensity (ξ) keeps contracting or relaxing until it reaches an extremum to switch memory value. Thus, the model is mechanistic and driven by the mechanics of the problem and it does not rely on any time-defined trigger to change state, it is based on maximum and minimum stretch of a lymphangion. A contraction cycle is defined by a contraction phase and a relaxation phase. If one considers that k c ω and k r ω are defined constants by equation (3.4), then the time constant can be written as the inverse of these products. In reality, they are not constant, since ω depends on the diameter, as shown in equation (3.5), so the value of these constants vary in time. Using figure 3.6, one can estimate the time spent between the upper and lower bounds of the contraction intensity. Furthermore, an example of a few relaxation-contraction cycles are provided by figure 3.13. If the terms k c ω and k r ω are regarded as constant, one can solve the differential equation (cf. equation (3.4)) and obtain that the time for relaxation is: t r = ecen-1 ecen+1 and time for contraction is: t c = ecen ecen+1 . By adding these two times, one obtains the period of the full cycle and dividing them, one gets the ratio of contraction-relaxation of a lymphangion.

Constitutive equation

The point of this section is to detail the relation between the pressure of the lymph and the diameter of a lymphangion. There are two types of wall forces here, one without muscle force and another one with muscle contraction. The first will be called passive pressure, the second one will be called contraction or active pressure. As listed in equation 3.1 (ii).

Passive pressure

The elasticity of the lymphangion walls leads to the existence of a passive component to the pressure, which is related to the diameter by the following non-linear relationship (3.9). This derivation is fully explained in the article of [START_REF] Bertram | An improved model of an actively contracting lymphatic vessel composed of several lymphangions: pumping characteristics[END_REF]:

P p = P d 4 15 12 exp D D d -1 -11 - D d D 3 (3.9)
This model only depends on the diameter. The constants are summed up in the table 3.2, furthermore they can be seen as gauges of the equation. Nonetheless, compared to the full curve shown in the article of [START_REF] Bertram | An improved model of an actively contracting lymphatic vessel composed of several lymphangions: pumping characteristics[END_REF], only the exponential term is used here as well as the cubic inverse term.

• The exponential term reflect on the elastin-collagen for high D Dd stretch, • The cubic inverse term is a collapse function at low pressure.

The figure 3.7 represents the passive pressure of the lymphangion as a function of the diameter as described in equation (3.9). This passive pressure prevents the lymphangion from dilating to a rupture point. It means that after a specific threshold related to D d , the diameter increases the pressure. Nevertheless, in the other case, if the diameter is small, the passive pressure should not have an important contribution to the constitutive law as shown in figure 3.7.

Active pressure

The active pressure is governed by equation (3.10), which relates it to three variables: the contraction intensity ξ, the thickness of the lymphangion wall τ and the diameter D.

P a = 2ξm 0 τ D -cD d cD d D . (3.10)
Where c is a constant that can be linked to the standard deviation of the stress-free diameter of the lymphangion. The thickness of the lymphangion τ is defined as :

τ = D -D 2 - 9 25 D 2 d . (3.11)
The contraction intensity and the thickness are only factors of the lymphangion diameter. One can observe that the active pressure can reach values three times higher than the passive pressure, using figure 3.8. Here, c represents a constant parameter related to the free stress diameter position, and also refers to as relaxed position of the lymphangion. It is obtained in terms of the D d constant.

Then by adding the contribution of the passive and active pressure get the figure 3.9.

Conservation equation

For this lymphangion modeling, the mass balance equation is used, and it is written like so:

∂V ∂t -Q net = 0. (3.12)
The partial derivative term represents the derivative of the volume V over time.

Q net shows the net flux entering the lymphangion. The volume V is computed in using the cylinder approximation, so one can write:

V = π D 2 [i] 4 (L -L valve ). (3.13)
With L the length of the lymphangion and L valve the length of the valve. Now that the reader is comfortable with the physics of this model and the approximation used here, we will now learn about the numeric aspect. The pressure here is measured in Barye and the diameter in centimeter. The range of diameter has been selected from minimum to maximum of these simulations.

Numerical implementation of a lymphatic network

Flow resistance

This model focus on specific diameter and pressure points, described figure 3.1. The physical approach used to the shallow water equations (Equations de Saint-Venant) like defined in this article [START_REF] Carson | An implicit solver for 1d arterial network models[END_REF]. We will consider three points for pressure and two points for the diameter. Concerning the pressure, P represents the pressure of the container, P v is the pressure just at the exit of the valve and P n represents the pressure at the entry of the valve. Since a zero-dimension model is used here, it is possible to represent a lymphangion using the electrical analogy, as in figure 3.10. Here, the spatial

P n[i-1] P v[i] P n[i] P v[i+1] P [i] P ext Q net Q 1 Q 2 Q 3 Figure 3
.10: Zero-dimension approximation of a lymphangion using the electrical analogy.

Dashed lines represent the limit of a lymphangion, diodes are the valve and capacitor represent the container/vessel.

discretisation is defined (using the subscript [i] ) based on the pressure drops between the elements of figure 3.10. If an equation is defined at a specific discretised time n and written like so [i,n] the subscript is implicit and is not written. It will only be written when different times are needed in an equation like so [i,n+1] , for example. Let us consider the point P v [i] , with K f [i] the conductance of the resistor and K v the conductance of the diode. We can write the following equations using the Hagen-Poiseuille law:

Q 1 = -K v[i] (P v[i] -P n[i-1] ), Q 2 = -K f [i] (P [i] -P v[i]
) and

Q 3 = -K f [i] (P n[i] -P [i]
).

(3.14)

K f is defined by:

K f [i] = 0.5 π(D [i] ) 4 128(L -L valve )µ .
It takes two K f to cross the entire lymphangion, hence the 0.5 coefficient in front of the fraction. K v is defined by:

K v[i] = π(D v[i] ) 4 128L valve µ with D v[i] = (1 -ζ [i] )(D [i-1] + D [i] ) 2 .
The variable ζ [i] is the valve state and is defined by the equation:

dζ [i] dt = (1 -ζ)k close (δP ) [i] if (δP ) [i] > 0, ζk open (δP ) [i] if (δP ) [i] ≤ 0.
Here, (δP ) [i] , is defined as

P v[i] -P n[i-1] .
Using the conservation of the flow, and since Q 1 is going to P v[i] and Q 2 is leaving it, one can add all the contributions to the flow and write:

Q 1 -Q 2 = 0,
Which leads to:

-K v[i] P n[i-1] + (K v[i] + K f [i] )P v[i] -K f [i] P [i] = 0. (3.15)
By applying the same discretisation to the point P [i] one obtains:

-K f [i] P n[i] + 2K f [i] P [i] -K f [i] P v[i] + Q net[i] = 0. (3.16)
Idem for the point P n[i] :

-K f [i] P [i] + (K v[i+1] + K f [i] )P n[i] -K v[i+1] P v[i+1] = 0.
(3.17)

Mass conservation

The mass equation is written for the lymphangion itself. It is written as follows:

∂V ∂t -Q net = 0 (3.18)
Q net represents the net flux entering the lymphangion, its value is only used in the matrix itself. The volume in this equation is estimated by the following decomposition ∂V ∂D ∂D ∂t . One can use the following approximation:

∂V ∂D = πL D 2 , and ∂D ∂t = 1 2δt (-3D [i,n+1] + 4D [i,n] -D [i,n-1] ). (3.19)
Here, the subscript ( [i] ) is for spatial discretisation and ( [i,n] ) is for the time. The second order upwind scheme is used for the time derivative and an analytical derivative for the diameter. This equation actually neglect the volume of fluid in the valve compartment.

Constitutive law

For this model, the constitutive law is used linking the pressure and the diameter but for different states of contraction. One part of this equation describes the pressure for a lymphangion with no muscle activity (relax state). It relates to the combination between the transmural pressure P p and the diameter of a passive lymphangion D. This part of the constitutive law can be found in this article [START_REF] Bertram | An improved model of an actively contracting lymphatic vessel composed of several lymphangions: pumping characteristics[END_REF].

P p[i] = P d 4 15 12 exp D [i] D d -1 -11 - D d D [i] 3 (3.20)
Here constants are P d and D d , representing the collapse pressure and the collapse diameter. These constants are defined in the parameter files. The second part of the equation focuses on the muscle activity (contraction intensity). To define the pressure due to the contraction, we need the wall thickness τ [i] of a lymphangion, the hoop strain σ [i] and the contraction intensity ξ [i] . The contraction intensity is defined by the following equation.

dξ [i] dt =      (1 -ξ [i] )k c ω if e min m 0 > ξ [i] , -ξ [i] k r ω if e max m 0 ≤ ξ [i] .
k c , k r , e min and e max are constants detailed in the previous section. ω is the frequency and is defined based on the diameter of the lymphangion D [i] . Then, one can derive the wall thickness τ such as:

τ [i] = D [i] -D 2 [i] - 9 25 D 2 d (3.21)
Here, m 0 is the wall stiffness, and the contraction pressure is defined by the following relation, using a linear parameter c to adjust the stress-free point:

P a[i] = 2σ [i] τ [i] D [i] = 2ξ [i] m 0 τ [i] D [i] -cD d cD d D [i] . (3.22)
These pressures are then introduced to the lymphangion at the point P [i] using the equation:

P [i] = P ext + (P a[i] + P p[i] ) + (D [i,n+1] -D [i,n] ) ∂(P a[i] + P p[i] ) ∂D [i] . (3.23)
The last equation used in this model is the conservation of mass, defined by equation (3.18), with the following discretisation as in the article [START_REF] Bertram | Incorporating measured valve properties into a numerical model of a lymphatic vessel[END_REF],

1 2δt (-3D [i,n+1] + 4D [i,n] -D [i,n-1] ) πL D [i,n] 2 -Q net[i] = 0. (3.24)
The five main equations are (3.15), (3.16), (3.17), (3.23) & (3.24). Which, can be summed up in the following system.

                         -K v[i] P n[i-1] + (K v[i] + K f [i] )P v[i] -K f [i] P [i] = 0 -K f [i] P n[i] + 2K f [i] P [i] -K f [i] P v[i] + Q net[i] = 0 -K f [i] P [i] + (K v[i+1] + K f [i] )P n[i] -K v[i+1] P v[i+1] = 0 P [i] = P ext + (P a[i] + P p[i] ) + (D [i,n+1] -D [i,n] ) ∂(P a[i] + P p[i] ) ∂D [i] 1 2δt (-3D [i,n+1] + 4D [i,n] -D [i,n-1] ) πL D [i,n] 2 -Q net[i,n+1] = 0
(3.25) Now that we have been through the main equations and their discretisation, we will see how to solve such system.

The system matrix

Using the system defined by equation (3.25), and with the following vector of unknown:

P n P D Q net P v T [i,n+1] , (3.26)
it is possible to write a system like AX = B, as mentioned in the previous section. In order to solve this system, we need to rewrite (3.23) and (3.24) to fit with this vector of unknown. So the equation (3.23) becomes:

P [i,n+1] + K a[i] D [i,n+1] = F a[i] , (3.27) with K a[i] = - ∂(P a[i] + P p[i] ) ∂D [i] ,
and

F a[i] = P ext[i] + P a[i] + P p[i] + D [n,i] ∂(P a[i] + P p[i] ) ∂D [i]
.

And the equation (3.24) becomes:

K m[i] D [i,n+1] + Q net[i] = F m[i] , (3.28) with K m[i] = 3 4δt πLD [i,n] ,
and

F m[i] = πLD [i,n] 4D [i,n] -D [i,n-1] 4δt .
The entire system is the equation (3.29) below:

      K v + K f -K f 0 0 -K v -K f 2K f 0 1 -K f 0 0 K m 1 0 0 1 K c 0 0 -K v -K f 0 0 K v + K f       [i,n]       P n P D Q net P v       [i,n+1] =       0 0 F m F a 0       [i,n]
(3.29) Once the matrix is assembled, it is solved. Then we replace the old variable by the new ones, and it is possible to write these data in a file or on the terminal stream. The matrix written in the equation (3.29), is written to be as symmetrical as possible.

Boundary and initial conditions

They are implemented in the time loop, as part of the right-hand side of the equation (3.16) only for the first lymphangion and the last of the chain. They are the only ones to have an unbalanced pressure drop equation. So one can write these equations:

K f [1] (-P v[1] + 2P [1] -P in ) + Q net = 0 K f [l] (-P v[l] + 2P [l] -P out ) + Q net = 0 (3.30)
Here, the subscript [1] means the first lymphangion of the chain and [l] the last one of the chain. Also, P in & P out refer to the selected pressure for the input and the output respectively. With the default inputs, one can find that half a minute, physics time, is enough.

To conclude this section, we just saw how the code works and how it solves the entire matrix system.

Integral quantities

Diameter average

In order to obtain the average diameter, a time integration is used here. Unless precised otherwise, it is counted from time zero and to the end of the physical time of the simulation, such as described by equation:

D = T 0 D(t)dt (3.31)
With D the diameter average, D(t) the average of all lymphangion diameter at time t. Furthermore, it can be compared to different specific diameters, such as the stress-free diameter.

Pressure average

Similarly to the diameter average metric, the pressure average is defined by a time integration. It is computed by the following equation:

P = T 0 P (t) dt (3.32)
In the equation (3.32) P , the pressure average, P (t) for a given time it represents the pressure average of all lymphangions, computed at its center.

Work estimation

Here, to compute the work properly we use the following approach:

W = P dV (3.33)
As lymphangions are regarded as cylinder, V = lπr 2 . P is the pressure of the lymphangion. Since we only compute the diameter of the lymphangion, and it depends in time, one gets the following simplification:

dV = 1 2 πlD dD dt dt (3.34)
However, there are two phases in time for a lymphangion, one spent contracting which contributes positively during time t c , and one relaxing contributing negatively during time t r . So we obtain the following relation:

W = tr 1 2 πlDP dD dt dt - tc 1 2 πlDP dD dt dt (3.35)
The term dD dt is easily computed by a simple difference explicit. t c and t r are defined using the contraction intensity variable shown in the equation (3.4).

Results

This section presents results from simulations run on a desktop computer1 . As well as Swansea University cluster called "Sunbird" via the supercomputing Wales portal2 . The code used here was implemented in Fortran 90 using four different files. All files holding functions related to:

• the mathematical part,

• the physical part,

• the input-output,

• the network aspect.

As well as, a main one holding the time loop and calling function form the other files. The matrix solver is using the MUMPS library [START_REF] Padua | Mumps[END_REF]. Furthermore, paralleling of independent simulation was obtained via the program parallel from [START_REF] Tange | Gnu parallel 20150322[END_REF].

Analysis of the behaviour of a simple vessel

We first analyse the behaviour of a single channel formed by a chain of three lymphangions such as represented in figure 3.11. This simple test case is a good set up to understand how the various models and approximations used can impact a lymphangion. Furthermore, this geometry is used on the first hand to look at specific variables such as the opening of the valve, the contraction intensity, as well as the mass flow and conservation. On the second hand, it is used to illustrate the behaviour of the model and the precision of the code. A simulation is run for 30 seconds of physical time with a time step of 10 -3 s. Imposed by the geometry shown in figure 3.11, there are two working valves, the first one is constrained by the input boundary condition. The boundary conditions were set to P input = 2 940 Ba or 3 cm H 2 O and P output = 4 900 Ba or 5 cm H 2 O. The initial conditions in pressure inside the vessel are identical to the input pressure.

Valve analysis

Here, a description is provided on how the model behaves in different regimes. Figure 3.12 shows the response of the valve in the central lymphangion. Using • a closing regime

• an opening regime

A third regime is implemented when the pressure difference is equal to zero, however such occurrence rarely happens. From a physical point of view, the valve closes under the action of the reverse flow velocity. But the two unknowns used in these simulations are the pressure and the diameter. So the pressure difference (δP ) is used. In this valve model, flow reversal is also assumed to be very quick when (δP ) changes sign so the associated relaxation constant k o and contraction constant k c are large. The value of this constant allows the switching time scale to be considered small from equation (3.2) defined previously. In figure 3.12 (iv) it is possible to compare the orifice diameter to a step function. Also, while the valve state is open, one can observe that the orifice diameter keeps varying, this phenomenon is linked to the equation (3.3) where this diameter is linked to the average diameter of the lymphangions surrounding the valve. No test has been conducted to modify the constant k o and k c to evaluate their range in order to obtain inertia type of behaviour. Figure 3.13: Evolution of the contraction of the second lymphangion in time. The graphs from top to bottom represent the contraction memory (i), the contraction intensity (ii), the pressure from contraction (iii) and the overall pressure (iv). Here, the contraction memory change from 0 (contracting) to 1 (relaxing). This allows the contraction memory to switch from contraction phase to relaxation phase, which then influences the contraction pressure.

The lymphangion pressure is shown to compare the contribution from the contraction pressure term.

Contraction analysis

Using the same simulation as in the valve analysis, one can plot the main variables for the contraction model in figure 3.13. To show the process and logic of this model, only the last ten seconds are displayed. The pulsating loop defined by exponential curves in the differential equation (3.4) is the variable (ii). Equation (3.4) shows that the state of contraction is a bounded function, and oscillates with values in the following interval: [9.9 × 10 -6 ; 6.0 × 10 -1 ]. These numerical value are set at the beginning of the simulation by emin m0 and emax m0 .

In a similar manner to the valve model with ζ 2 of figure 3.12, here γ 2 , is the gate function, with zero being the contracted state and one being the relaxed state. From a physics viewpoint, if one compares the curve of the contraction memory and the ones from the valve state, the timescale is approximately the same as there is three cycles in both models. Since the contraction is forced and imposes the pressure, the valve cycles will synchronise with the contraction cycle after a while. Here, for example the period is of approximately 2 s.

While this contraction model can be quite periodic, the previous model for the valve is not. If they are synchronised then it suggests that all the lymphangions are coordinated since the valve can be seen as the interface between two following lymphangions. However, this is not always the case, it has been commented in the following section. This non-coordination of contraction is explored for longer vessel later in the thesis.

Time scales of contraction and relaxation are defined by constants k c ω and k r ω. It is important to notice, ω depends on D, this relation makes the model non-linear. This coupling with the diameter D causes variation on the time scale of the contraction cycle. Moreover, trends can be observed by varying parameters such as boundary condition pressures and initial conditions. Again with the simulation used in the previous paragraph, it is possible to observe the contribution of the constitutive equation of this model in figure 3.14. For the contraction phase, it is defined as the sum of the passive and active terms defined in equations (3.9) & (3.10). Furthermore, we compare it with the total pressure of the lymphangion and the diameter of lymphangion.

Constitutive equation analysis

The main contribution on the pressure is achieved via the constitutive equation.

Figure 3.14 (iii) represents the volume vs pressure inside the second lymphangion. It shows almost three cycles, of contraction and relaxation. The first one starts at the point A then goes to B and finishes at C. A represents the start of the contraction, it has the largest volume before the contraction. B shows the start of the relaxation, an isovolumetric pressure drop occurs then the volume increase at the same rate until the next contraction, from the point C. The second cycle is completed between points C, D and E. A third cycle can be observed as well between E, F and a point situated between C and E. The volume span is shorter for the second cycle compared to the first one. Nonetheless, the work produced in these two cycles is approximately the same, since the second cycle provides a stronger contraction with higher peak in pressure. In order to test the model against reality, it is important to obtain experimental data. As this aspect is not covered in this thesis, we will just use data gathered by [START_REF] Michael | Intrinsic increase in lymphangion muscle contractility in response to elevated afterload[END_REF] & [START_REF] Bertram | Consequences of intravascular lymphatic valve properties: a study 135 of contraction timing in a multi-lymphangion model[END_REF]. However, the units used in those papers are not always CGS unit system. Data used here was digitised by engauge digitizer3 . They used laboratory rats to conduct their experiments. By extracting and cropping the lymphatic vessel to a two-valve segment as shown in figure 3.15, they inserted glass micropipette at each end of the vessel. This set up allows to measure the input and output pressure via a low-pressure transducer. To obtain the diameter values, they used a monochrome camera wired to a computer tracking the wall movements.

Comparison with experimental results

P out P in P L central lymphangion
Nevertheless, two experimental devices were designed to simulate contraction from the wall and from the valve. More details are provided in Appendix A. In this section, different results obtained in experiments are analysed.

Parameters tuning

First, the relevant physiological ranges need to be defined in order to vary parameters properly. The table 3 e max e min D d P d m 0 α c 1 c 2 f min 3 × 10 5 10 -8 2.9 × 10 -3 36 10 8 1.9 1.9 2.5 × 10 -1 10 -2 2 × 10 4 10 -4 2.9 × 10 -3 35 10 8 2.0 2.4 2 × 10 -1 5 × 10 -2 relaxation phases, as it is a ratio, it is therefore dimensionless. Table 3.4 lists all the parameters used to tune the model to experimental data gathered on different articles. These were manually tuned using batches of simulations to scan the parameter space.

Comparison between simulation and digitized data

Data set 1 This data set was gathered from figure 1 of the article Davis et al. [2012]. After running batches of simulations and varying the parameters listed in table 3.4, a behaviour comparable to experimental results in terms of frequencies and amplitudes was obtained with the parameter values listed in table 3.5.

Furthermore, in order to model the boundary conditions properly from 1 cm H 2 O to 6.5 cm H 2 O then drops quickly back to its original level. The growth slope of used in the simulation is estimated as 2.75 cmH 2 O/min.

Pressure is better captured than the diameter. Difference may come from the Figure 3.17: Diameter (cm) of second lymphangion vs time (s) for the simulation (i) and the experiment (ii). In terms of frequency, one can observe faster contractions when the pressure difference is larger and amplitude is reduced. Maximum of peak values are constant in the simulation while they decrease in the experimental data.

P out P in 1 2 3 (i) (i) t(
constitutive equation or the contraction cycle. This difference might also come from the irregular pulsation cycle of the experimental data, while the simulation imposes more regular pulsations. The maximum diameter from the simulation remains constant while the experimental data appears to have two plateaus, one before and one after 250 s. The first one is close the 9.5 × 10 -3 cm, and the second one with faster contraction is lower than the initial one: 9 × 10 -3 cm. Furthermore, regarding the shape of the diameter evolution for one contraction cycle is not properly computed. Indeed, the simulated diameter overvalues the time spent for each local minimum before 250 s in figure 3.17.

With regard to the frequency, figure 3.18 shows that simulation data have higher frequency than experimental one. As we saw before, the experimental data hold two different regimes, one before 250 seconds and one after. For the first one, experimental data gives approximately 7.2 contraction cycles per and the digitized data (ii). Here the envelop matches well between the two graphs. The frequency trend is properly captured as well. minute (ccpm), while the simulation is closer to 11.4. On the second phase and until the output, boundary condition gets back to initial level, from 250 s to 330 s, there is approximately 14.25 ccpm in the experimental data, compared to 19.5 for the simulation. In both cases, the number of contraction per minute is overestimated in the simulation. Even if this difference is large, what matters here is to capture the envelop of the pressure oscillation. Data set 2 Here, we focus on experimental results provided in figure 9 f min 2 × 10 4 10 -4 2.9 × 10 -3 35 10 8 2.0 2.4 2 × 10 -1 5 × 10 -2 3 × 10 5 10 -8 2.9 × 10 -3 36 10 8 1.9 1.9 2.5 × 10 -1 10 -2 comparing the two sets of data, few parameters need to be changed, the one which holds most variation is e min . In order to model the boundary conditions properly, we follow the input and output described in the paper. In this data set, the pressure values inside each lymphangions are not measured by article's authors, so we only focus on the diameter. However, it is possible to observe a difference in the growth starting time. For the simulation, it happens straight after the change in boundary conditions, so approximately around 210 seconds. For the experimental data, there is a 30 seconds delay, as the diameter and the frequency switch.

Regarding the frequency, there are multiple observations:

• For simulation data the average is 10 ccpm for all lymphangions except the last (5th) one.

• The last lymphangion of the simulation data has 18 ccpm before the ramp up but goes back to average when the boundary condition pressure falls back to 2 cm H 2 O of pressure difference.

• Experimental data before the variation of the boundary conditions:

-Lymphangions are on average 10 ccpm.

-The two initial lymphangions appear to contract less.

• Experimental data after the variation of the boundary conditions (t > 440 s):

-Lymphangions are on average 4 ccpm.

-The two initial lymphangions have half the frequency of the third lymphangion.

-The last two lymphangions only contract once during a minute.

-Only the third lymphangion recover approximately its previous frequency.

In the end, the simulation tend to overestimate the frequency after 440 seconds of simulation. Furthermore, before 210 seconds the experimental data have non-stable contraction behaviour, and such behaviour cannot be produced by simulation.

In terms of diameter, the simulations rise up to 0.015 cm during peak pressure on all lymphangions. The same behaviour is observed in experimental data, except the peak diameter is slightly lower and around 0.013 cm. The third lymphangion has the largest amplitude in the experimental data while it is the first one for simulations. Furthermore, diameter is similar for all lymphangions in the simulation while figure 3.4.2 shows that first two lymphangions have small amplitude compare to lymphangion three to five in experiments.

To sum up, there are good agreements between the model defined here and the experimental data in terms of order of magnitude. But the frequency in the experimental data is too chaotic to be captured properly. Numerical data tends to overestimate the number of contraction. Furthermore, mouse lymphatics network is analysed and fit the model according to the experimental data. As we try to model a full human leg network, the similarities between such rodent and human on their lymph system is not well studied, some literature exist on this topic and the length and time scales appear to be the main factors of difference [START_REF] Mestas | Of mice and not men: Differences between mouse and human immunology[END_REF].

Three-vessel networks

Here we introduce bifurcation and confluence into the vessels. This type of network has partially been studied by [START_REF] Jamalian | Network scale modeling of lymph transport and its effective pumping parameters[END_REF]. The secondary lymphatic network found in a leg skin is ramified and shows various bifurcations and confluences. In fact, all the capillaries merge into initial secondary vessels and their number are greater than lymph node total input. The overall secondary vessels network has a converging structure, but still displays some bifurcations, in this type of network with many lymph inputs and few outputs. In order to understand the dynamics of the network, we run different simple geometries such as a confluence and a bifurcation. In the following geometries, e max e min different number of lymphangions per vessels are considered. Moreover, a new variable is added for the boundary condition: δP BC = P out -P in . Both input and output have steady pressure value during the simulation, unlike previous cases. For all simulations, the input boundary condition is defined as 5 880 Ba (6 cm H 2 O). First, we are going to study the confluence as shown on figure 3.22 on the left.

D d P d m 0 α c 1 c 2 f min 2 × 10 4 10 -4 2.9 × 10 -3 35 10 8 2.0 2.4 2 × 10 -1 5 × 10 -2
The metrics used to quantify the behaviour of the system are the average work per contraction and per lymphangion, as well as the number of contraction.

Here, the parameters used for the contraction model of this simulation are summed up in table 3.7. Furthermore, all the simulations are ten minutes long. A minute is used at the beginning to stabilise the system and is not taken into account.

The key parameters are the pressure difference on the imposed on the boundary condition (δP BC ), and the number of lymphangion per vessel (n). It is important to notice that whatever the number of lymphangions in a vessel, their length does not vary, and is defined as 10 -1 cm, and that all vessels of the network get the same number of lymphangions. For example, if one considers a vessel with 11 lymphangions then the total vessel length would be 1.1 cm. By varying n it will allow us to understand under which condition the system behaves more efficiently. Main interest is to observe the influence of pressure difference over the boundary conditions. The main two metrics used here are average work per lymphangion and average contraction per minute. To estimate the work metric, equation (3.35) is used, then the result is averaged over total time and lymphangions per vessel. Same process is used for the average contraction per minute.

Confluence

In this study, we will consider the number of lymphangions per vessel as a variable. However, it is the same for all the three vessels of the network. It is important to note that only two vessels graphs are shown (afferent top & efferent), since the geometry has a top-down symmetry. Furthermore, the boundary conditions are symmetric as well, so the curves are almost identical, between the vessel v c1 and v c2 . The subtle modification between the two afferent vessels will not be considered here. 3.23 (i) shows that the afferent vessels tend does not balance the work between each lymphangion properly. In fact, the lymphangions closer to the merging node produces more work than the one close to the input of the vessel. It is less the case in the (ii) graph, as one can observe that the work is approximately the same for all lymphangions from the confluence node to the output. .24 (i) can be distinguished into two main groups: vessels with great number of lymphangions (n ≥ 7) and vessels with less than seven (n < 7). The first one grows more or less linearly with the pressure difference on the boundary conditions. However, the second group behaves differently. Indeed, for most cases when δP BC is set to 4 cm H 2 O, a maximum appears and configuration with three lymphangions displays a peak work. Also four lymphangions configuration peak work is shown for δP BC of 7 cm H 2 O.

In figure 3.24 (ii), the overall behaviour of this vessel is quite different from the previous one. It shows that whatever the number of lymphangion per vessel, the greater δP BC the larger work is provided in order to push lymph. It is possible to observe as well that for large δP BC (7 cm H 2 O and above) the behaviour is quite similar. A simple explanation could be that more lymphangions per vessel result in lower to individual pumping work. Besides, for large value of n and small values of δP BC , the average work appears to be constant, which is the base level the system remains when small difference of pressure are applied, such as observed in the previous section. lymphangion. It allows to observe whether the lymphangion would have a peak in activity, and draw conclusion whether the work load is spread in all vessels of if afferent or efferent vessels are more solicited. We observe that whatever the number of lymphangion in a vessel (n), averaged frequency grows with the pressure difference δP BC for vessel v c3 . This is also observed for the vessel v c1 , and one can notice that the frequencey does not vary with the lymphangion numbers (n). The less lymphangion in a vessel, the more often they contract for afferent vessel to the confluence. Furthermore, vessel v c1 has an overall average contraction frequency higher than v c3 . These differences between v c1 and v c3 lead us to the conclusion that efferent vessels in a simple confluent network are more active than efferent ones, whatever the size or the pressure difference on the boundary conditions. The word active here relates to an average frequency, and a higher average work output.

Bifurcation

In this section, the main focus is the bifurcation network which is represented in figure 3.22 (ii). We reproduce the same geometry set as in the previous section except with a bifurcation network. Similarly, when the number of lymphangion per vessel is set to n, all vessels are made of the same n number of lymphangion. In figure 3.28, very similar results are observed compared to figure 3.24, but the scale is reduced by half. It appears the lymphangions are less active than in the confluence type of network. However, in this configuration there are two out vessels while on the confluence network, there is only one single out vessels. One can notice that the more lymphangion in the vessel the less an individual lymphangion produces work. Also the greater the boundary condition δP BC the more lymphangions will produce work.

Figure 3.29 (i), shows even more clearly than in figure 3.30 that vessels with more than seven lymphangions are more efficient with pressure gradient near 7 cm H 2 O. Moreover, three to five lymphangions vessels are more efficient with the lowest pressure gradient. For high pressure gradient all vessels appears to show linear pattern. As a matter of fact for high δP BC , for 3 cm H 2 O the average work per lymphangion multiplied by n grows by 24 × 10 -3 erg, so a slope of 1.25 × 10 3 cm H 2 O erg -1 approximately for all vessels.

For the efferent vessels of figure 3.29, the average work per lymphangion multiplied by n grows for δP BC larger than four. Furthermore, one can fit the numerical results against a mathematical function. Two types of functions were The best fit for the power law provides n ≈ 3, so the average work per lymphangion multiplied by n grow by a third of the boundary conditions pressure difference. But if one consider figure 3.28 without multiplying by n then for the same high δP BC and high range of n, the slope is more linear, while for low range of n a power law or an exponential is more suited. In figure 3.30, the tendencies are quite different from the ones observed in the previous network. Here, the amount of work appears to be decorrelated to the pressure difference of the boundary condition. But the same trend is observed as the number of lymphangions per vessel is inversely proportional to the work given.

The average contraction period for vessel (v b2 ), is larger than its equivalent in the previously considered network. This metric also shows that the lymphangions are contracting more and thus more active. To put this in perspective with the work provided, it appears that lymphangions contract more, but the average work is not as great as for the confluence network.

To conclude, the bifurcation network tends to provide better balance in the work provided by both afferent and efferent vessels. Nonetheless, the efferent vessels have a higher frequency, and this can be translated into a saturation and difficulty for the vessel to meet the boundary condition. The number of lymphangion per afferent vessel does not modify its behaviour, and only impacts its frequency for pressure difference lower than 7 cm H 2 O.

Leg skin networks

Introduction In this section, the dynamics of the model in a realistic lymphatic network geometry is studied. This geometry was acquired from original drawings of a leg skin lymphatic network by Sappey [1874]. Initially, we describe how the digitisation process works, then a reflection on the general structure of the network is proposed. Then, we inspect different vessel regime and contraction wave of different lymphangion vessels. 10 6 10 -8 10 -3 35 10 8 1.9 1.9 2.5 × 10 -1 5 × 10 -2 input and output are defined by selecting the area used for the simulation. As detailed in the introduction, the lymph flows from the interstitial fluid back to the blood circulation near the heart. So if we apply this into the original images shown in figures 3.31 & 3.32, the lymph flows approximately from left to right, as it enters the network from the foot or capillaries and exits the network by accessing to lymph nodes of the groin area. As we study the leg network, the simulation domain is delimited by the red lines at the ankle and the red line in upper part of the leg, near the groin area. Then, bifurcations, confluences, inputs and outputs are identified and labeled. In the network of figure 3.31, there are 72 inputs for 25 outputs, so approximately three times more inputs.

In figure 3.32 one can observe 46 inputs and 26 outputs, so the ratio is closer to two times. As a consequence of these ratios, networks are dominated by confluences, which reflects the collecting function of the network. There is still a small, but significant number of bifurcations though as seen in table 3.9. In order to convert the network from figure for its simulation, it is necessary to identify the beginning and the end of each vessel. So the limit of the context is defined in both figures 3.31 & 3.32 by the red line. The dimension chosen for the leg are 100 cm of length and 30 cm half the circumference. Furthermore, the bifurcation are referenced as blue dots, and confluence as yellow dots. The coordinates are used to then evaluate the distance from the beginning to the end of a specific vessel. Knowing the distance, the number of lymphangions in a vessel is then determined automatically by the average size of a lymphangion provided in the input file. A lymphangion is estimated to be 0.3 cm as define in table 3.3 as L. Such information allows to build the entire network. Since the position of lymphangion is not important for the computation it is not used in the code. Data set to build the networks for these simulations are provided in the appendix C.

In terms of boundary conditions, the inlets are set to a pressure of 2 940 Ba, and the outlets to 4 900 Ba. This applies a δP BC of 1 960 Ba across the whole network. The parameters used to model a lymphangion are all defined in table 3.8. However, it is used again for post-processing such as figure 3.37 & 3.38. It is interesting to note, position of vessels might as well be used to include the effect of gravity for a later integration. Chaotic behaviour Table 3.10 provides the position of the 61st vessel of the first drawing as well as its input and output. Figure 3.33 represents the 61st vessel which is an input vessel going from the middle top of the physical domain to a merge point. It is possible to notice the merging of 56th and the 61st into the 72nd vessel. As one can observe, between 20 s and 30 s, the last few lymphangions (11 to 14) do not contract and remain still, except the last one (15). These lymphangions (11 to 14) are still in the initial equilibrium state, which they exit as the lymph is carried out of the vessel. In this initial state, a lymphangion remain still and does not contract cyclically. For this vessel, no train of pulsation appears to synchronise all the lymphangions. The average number of contraction for a minute is about 15 to 20 per minutes which is relatively high compare to literature findings. For a human the thoracic duct the frequency of contraction is found to be in between 5 to almost 15 per minutes according to [START_REF] Telinius | Human thoracic duct in vitro: diameter-tension properties, spontaneous and evoked contractile activity[END_REF]. But since these lymphatics are trunks, and we are focusing in the secondary collecting network one can use the Womersley number which is in our case closer to 0.1 than unity, as precised in [START_REF] Moore | Lymphatic system flows[END_REF]. This dimensionless number represents the ratio between transient inertial force to the viscous force, and is defined like so:

Wo = L ω ν (3.36)
Since the diameter of a secondary lymphatic vessel is smaller in diameter than a trunk vessel, then the Womersley number should be smaller, which increase the pulsation of the transient inertial force term. Here, viscous force dominates over the transient inertia, while in trunks they balance each others. From figure 3.33, it is possible to compute ≈ 0.12 for the Womersley number. All in phase behaviour Table 3.11 gives position, input and output of the 17th vessel of first drawing as well as 103rd vessel of first drawing network. Figure 3.34 represents a vessel on the beginning of the network near the ankle from the first drawing. It is an input vessel like the previous case, here it leads to a merge node going to the 24th vessel. However, here one can observe that all the lymphangions are in phase with each other after 50 seconds. Before these 50 s, the group formed by the last three lymphangions are behaving chaotically. It is possible to follow the path of wave of contractions and see how it spreads in a vessel. The red line in figure 3.34 shows the propagation of the pressure pulse across 15 lymphangions in 10.5 seconds, which indicates a delay of 0.7 s and a propagation velocity of 0.28 cm s -1 Figure 3.35 display the same type of synchronised behaviour, while being on the other side of the leg close to the groin lymph nodes. As all the previous vessel analysed, this one ( 103) is an input vessel going to a merge node, except it is longer by three lymphangions. After fifty seconds it is possible to see three waves well spaced, but with slower speed. In red, figure 3.35, shows the waves propagating from the third lymphangion. Furthermore, I decided to not use the last lymphangions to estimate the velocity for the last wave as the sixteenth one contracted almost as soon as the fifteenth one. So the two waves here have a propagation velocity of ≈ 0.31 cm s -1 and ≈ 0.2 cm s -1 .

Appearance of synchronisation behaviour

In figure 3.36 represent the pressure for all lymphangion of the 72nd vessel, it results from the confluence of two vessels and lead to a bifurcation. Table 3.12 provides the position of the 72nd vessel of drawing 2 as well as its input and output. There are couples seconds, a wave appears from the 6th lymphangion, but the following ones are broken by the 15th lymphangion. Some similar events occur later with third lymphangion but the relaxation phase is separated by the ninth lymphangion. Couple of reasons make this vessel different from the other cases studied. On the first hand the number of lymphangions which is larger than previously shown vessels. On the other hand, its position in the network is different as it is not an input vessel, but it is fully integrated in the network. Overall view From figure 3.37, it is possible to observe screenshots as the simulation evolve. Every dots on figure represent a lymphangion. For some vessel, lymph never reach them and thus never start to contract automatically. The simulations did not run for enough time to observe a privilege route for lymph. However, few routes showed to be sub-optimal. This is not pushed further since the code build bifurcation such that it might lead to lymph spillage, and obviously need further investigation.

Concluding remarks

We have seen that this lymphangion model can compute the lymph dynamics of an entire secondary network. Different behaviour has been identified, but the equilibrium state has not been reached. It would be interesting to let the simulation run longer than 80 s. From these longer simulations, some privilege route should appear. Furthermore, inlet and outlet flow can be computed to compare with clinical data. Some could suggest a chaotic approach in order to deeply understand the dynamic of bifurcations Abstract -The work achieved here proposes a two-dimensional model of a lymphangion, and study the behaviour of lymphocytes in the lymph flow. First, by giving the fluid and solid equations of the simulation model used. How they are coupled together with all the technical details related to their implementations in MATLAB. Then, we study the behaviour of a group of lymphocytes in a lymphatic channel with valves using an FSI (Fluid-Structure Interaction) code. We define a lymphangion geometry and compare the displacement of lymphocites in different scenari, first with or without valve, pulsations and/or movable walls. The lymph flow in the secondary network is similar to blood pulsating flow.

Numerical method

In this section, we aim to provide a description of the equations, models and numerical methods used for the 2D simulations. We will first present the computation techniques for the fluid, followed by those for the solid. The fluid-structure interaction scheme is given after this description of the computation techniques. A lymphangion as defined in figure 1.1 in the introduction contains fluid and cells. It carries lymph which for the purpose of this thesis will be considered as a classic incompressible Newtonian fluid. However, this hypothesis is not always true as the density of biological cells varies and modifies the fluid viscosity and possibly its visco-elastic properties. Another known alteration of rheology in the lymph network occurs in the small intestine during digestion. This altered lymph is called chyle and contains lipids which modify the viscosity Bragdon [1958]. Incompressibility is a property of the flow that allows simplification in the fluid solver such as the density is considered as a constant. A newtonian fluid correlates shear stress (τ ) and local strain rate of the fluid, in such a way that they are proportional: τ = µ(∇u + ∇u T ). On the other hand, the valve lymphocytes and walls will be considered as incompressible hyperelastic solids following a Neo-Hookean model. Such model links the Cauchy stress tensor (σ), with the pressure (p) and the Cauchy-Green deformation tensor (B). So one can write: σ = -pI + 2C 1 B, with (C 1 ) a shear modulus constant. Furthermore, the density of the fluid will be written ρ.

FSI (Fluid-Structure Interaction) codes are multi-field problems since the fluid and solid share a moving interface. There are different strategies when it comes to choose the best support to capture the interface. A classic formulation such as ALE (Arbitrary Lagrangian Eulerian) uses a moving spatial domain to resolve the solid position and movement as described in [START_REF] Maury | Fluid-particle flow: a symmetric formulation[END_REF]. This fluid & solid domains requires more computational power as its boundary conditions and mesh need to be computed at every time step [START_REF] Bavo | Fluid-structure interaction simulation of prosthetic aortic valves: Comparison between immersed boundary and arbitrary lagrangian-eulerian techniques for the mesh representation[END_REF]. Also, this method is more suited for limited deformations of the solid geometry. Another formulation called IBM (Immersed Boundary Method) developed for both incompressible fluid and solid shows better results for high deformation of structure, which is important in biomedical applications. In fact, it uses a different way to share information between the fluid and the solid, such that the solid interpolates the fluid solution in order to compute its deformation, displacement and forces. The forces are then added back to the Navier-Stokes equations to include the effect of the solid. However, this interpolation can lead to spreading and excessive diffusion of the numerical solution. But in this code the solid is treated as a potential energy functional completely immersed. Also, different solid computational point1 functions has been considered in order to reduce the numerical diffusion.

The FSI method used here is introduced in the article by [START_REF] Gil | The immersed structural potential method for haemodynamic applications[END_REF] and is called ISPM (Immersed Structural Potential Method) and is a variation of the IBM. The flow chart shown in figure 4.1 shows the organisation and links between the main steps in order to run the code. The starting point the chart is In green rhombi are the conditions occurring in the code logic. Here, t represents time, ∆t is the time-step, Err(fFSI) =

|f n FSI -f n-1 FSI | |f n FSI |
the error function and k tolerance is the tolerance constant.

the definition of the simulation. It is a MATLAB file from which one can choose the methods for boundary conditions, size for the physical domain, position of the structure and other elements related to the geometry. Then the fluid is computed using a predictor-corrector method which consists in computing the viscous term then the inertia term which are then summed into a predicted fluid velocity variables. Next, the boundary conditions are applied on the velocity fields. Then, these new fluid variables are used to compute the pressure. The gradient of the pressure is used to update the predicted velocity into a corrected velocity.

After solving the fluid, the structure is computed, first by interpolating the velocity values of the fluid using kernel functions. From there, the velocity gradient is computed in order to get to the Kirchhoff stress tensor which will give the FSI forces to be applied to the Navier-Stokes equation. Then the sodid computational points are displaced based on the interpolated velocity. Using the FSI forces field, the relative error of the forces is compared to a tolerance constant at every iteration. A new time step is computed when the relative error is below a defined tolerance. Now that the key steps for the ISPM code have been introduced, the detail of the main equations is provided in the following subsection of the thesis.

Fluid computation

Mesh

The computational domain (Ω) is defined by a rectangular fluid mesh, also referred to as the structured grid. It is composed of uniform rectangular cells. This mesh holds the values of velocity along the x-axis and y-axis, as well as the pressure field. For stability reasons, the computational point of pressure and velocities are shifted as presented in figure 4.2. In this figure, the lines form the computational pressure cells. This method of staggering the pressure points and velocity is also called Marker And Cell (MAC). It was first proposed by [START_REF] Francis | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF]. These computational points are also defined for the u and for v, which will be useful for the Quadratic Upstream Interpolation for Convective Kinematics (QUICK) scheme.

u i,j u i+1,j v i,j v i,j+1 P i,j

Boundary Conditions

In order to apply Boundary Conditions (BC), an additional layer of fluid cells is added at the border of the simulation domain, these cells are called ghost cells. They are not taken into account in the results but play key role in the calculation of the solution. From a numerical point of view, they allow all the multiple differential operators to provide results for all the cells of the simulation domain. In our case they account for two cells in addition to the computational domain just for the BC. The number of additional cells depends on the scheme used by differential operators. Furthermore, the value of pressure and velocity held by these cells is defined before solving the physics thus they are known a priori and influence the results. There are various types of BC, for the velocity field in computational fluid dynamic. The main types of BC used in this thesis are:

• Dirichlet

• Neumann

In order to have Dirichlet type of BC, all velocity values in the BC cells are defined and fixed at the beginning of the simulation. So one can write:

u(x) = f D (x) with x ∈ ∂Ω.
With ∂Ω the boundary of the domain, and f D a function defining the velocity at the BC. A Neumann BC is different as it depends on the value of u near the border of the simulation domain, since it is function of a differential term in space. It can be written:

∂u ∂n (x) = f N (x) with x ∈ ∂Ω.
Here n is the normal to the boundary ∂Ω and f N the function describing the flux of velocity. Furthermore, a periodic BC is used in one of simulations presented in this thesis. More details are provided in the definition of the geometry in the specific section later in the document.

Fluid numerical method

In this section, is a description of the major steps of the numerical formulation as described in figure 4.1, in order to solve the fluid fields, which is the pressure, and the two components of the velocity. Also, the fluid velocity vector field is u = (u v) T .

Time progression

This code uses the predictor-corrector method, which are explained in [START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF]. Using the Navier-Stokes equation one can obtain for the predicted step:

u * -u n δt = - 3 2 (u n • ∇)u n + 1 2 (u n-1 • ∇)u n-1 + 1 ρ (µ∆u n + f n + f n ine ). (4.1)
The name of the time discretisation scheme is called Adam-Bashforth scheme.

Here, u n is the velocity at the n-time step, u * is the predicted velocity and f n are the structural forces applied to the fluid. Furthermore, f n ine represents the inertial forces which are only taken into account if the density of the fluid and the solid are different. At this step the boundary conditions for the fluid are imposed as:

u * | ∂Ω = u n+1 BC . (4.
2)

The correction step can be written as follows:

∆p n+1 = ρ δt ∇ • u * . (4.3)
Here p is the pressure field and ρ is the fluid density. The Poisson equation ( 4.3) is solved using a Fast Fourier Transform (FFT) solver in order to improve efficiency. The boundary conditions are then applied as:

n| ∂Ω • ∇p n+1 = 0. (4.4)
Finally the velocities are corrected:

u n+1 = u * - δt ρ ∇p n+1 . ( 4

.5)

Viscous discretisation The diffusive term of equation ( 4.1) is written as follows:

(∆u) i,j = µ ρ u i-1,j -2u i,j + u i+1,j (δx) 2 + u i,j-1 -2u i,j + u i,j+1 (δy) 2 . (4.6)
Speaking of computational efficiency, the mesh size (δx & δy) is constant as well as the viscosity (µ) and the density (ρ) so this operator can be vectorised, for i and j in the entire simulation domain, except the ghost cells. The equation which defines (∆v) i,j is similar and thus is not written here.

Inertia discretisation

The QUICK formulation requires to know the velocity normal to the computational cell, represented in figure 4.2. If one consider the point u i,j the orthogonal velocities of this cell are defined as follows:

ūi+ 1 2 ,j = 1 2 (u i,j + u i+1,j ), ūi-1 2 ,j = 1 2 (u i,j + u i-1,j ), ūi,j+ 1 2 = 1 2 (v i,j+1 + v i-1,j+1 ), ūi,j-1 2 = 1 2 (v i,j + v i-1,j ). (4.7)
ū represents the averaged velocity in between two neighboring computational velocity points. Knowing the velocities normal to the cell, one can apply the QUICK formulation to obtain the interpolated velocities û :

ûi+ 1 2 ,j = 1 8 (3u i+1,j + 6u i,j -u i-1,j ) if ūi+ 1 2 ,j > 0 1 8 (3u i,j + 6u i+1,j -u i+2,j ) else. (4.8) And, ûi-1 2 ,j = 1 8 (3u i,j + 6u i-1,j -u i-2,j ) if ūi-1 2 ,j > 0 1 8 (3u i-1,j + 6u i,j -u i+1,j ) else.
(4.9)

In the same manner, one can obtain the formulation for the y-axis, and find the following velocities, ûi,j-1 2 and ûi,j+ 1 2 . The convective term of the equation (4.1) is then written as follows in the code:

(u i,j • ∇)u i,j = ūi+ 1 2 ,j • ûi+ 1 2 ,j -ūi-1 2 ,j • ûi-1 2 ,j δx + ūi,j+ 1 2 • ûi,j+ 1 2 -ūi,j-1 2 • ûi,j-1 2 δy .
(4.10)

The same procedure is applied to compute the v i,j contribution to the convective term. For i and j in the entire simulation domain, except the ghost cells. In terms of vectorisation it is still manageable but the condition on the equation (4.8) and (4.9) required to create list of boolean in order to "unroll"2 the for-loop properly.

Solid computation

Points -Kernels

In a same way that the fluid needs a mesh, the solid structure need computational points to evaluate its position. These points are also referred to as kernels then a cloud or collection of kernels form a solid structure. The formulation for kernel is defined in [START_REF] Peskin | The immersed boundary method[END_REF], which is the same as the one used here follows these equations :

φ(r) =      1 8 (5 -2|r| --7 + 4(3 -|r|)|r|) if 1 ≤ |r| ≤ 2, 1 8 (3 -2|r| + 1 + 4(1 -|r|)|r|) if |r| < 1, 0 else. (4.11)
The φ function is shown in figure 4.3, first as a function of the radius r and then in two dimensions.

These kernels are the support of the solid structures, each one hold the position of the solid as well as its density They are usually defined as primitive shape such as a square or a circle at the start of a simulation. It is then possible to adapt their geometry by applying geometric transformation. In the remainder of this section, N will refer to the number of solid structure points, and x s the point coordinates of a kernel. Any variable with subscript ( s ) will refer to the collection of solid structure points defining a solid in the simulation. Besides, if subscript ( n ) is not written, it is implicit.

Boundary Conditions

In terms of boundary conditions for the solid, The solid structure depends on the fluid computational mesh. So the solid BC are defined by the underlying fluid conditions as detailed in figure 4.1 The solid kernels are used to compute the structure physical position, displacement and velocity through interpolation of the fluid mesh. If, one would like to get the solid attached to a position in the simulation domain, it would then need to define a patch of zero velocity for the fluid which will then interpolate into the solid.

Solid numerical method

In this subsection the structural displacement and forces computed by the algorithm from the interpolated fluid velocity field is detailed.

Interpolation A kernel approach is used to compute the interpolation of fluid velocity, as all the kernel use an interpolation function. In the article [START_REF] Gil | The immersed structural potential method for haemodynamic applications[END_REF], few different ones where tested. This approach defines a distance of support for the kernel, to extract velocity from the fluid mesh. Then the velocity of a kernel is defined by:

u s = Ax∈Ix u Ax φ Ax (x s ). ( 4

.12)

With A x being a fluid cell edge within the kernel support written I x . Furthermore, φ is the kernel function defined in the previous subsection. The same is achieved to obtain v s using A y .

Velocity gradient

The solid velocity gradient is computed as:

l s = [∇ s u s ∇ s v s ] T (4.13)
Numerically, the gradient is applied to the the kernel function φ from the equation (4.12).

Deformation gradient tensor Then, in order to obtain F s the deformation gradient tensor, we will use this formula:

∂F s ∂t = l s F s . (4.14)
Which lead to for the n + 1-th time step:

F n+1 s = e ∆tl n s F n s . (4.15)
To properly estimate the exponential term, one need to separate l s into its symmetric component d s and skew-symmetric w s . This lead to:

e ∆tls = e ∆tds e ∆tws . (4.16)

From a numerical point of view, the d s term is solved using spectral decomposition. And the w s term is solved by Hughes-Winget's update method as developed in Hughes and Winget [1980].

Deviatoric Kirchhoff stress tensor

The stress tensor is obtained using an incompressible hyperelastic Neo-Hookean material constitutive equation:

ψ s (C) = G 2 (I c -3) and I c = λ 2 1 + λ 2 2 + λ 2 3 = tr(J -2/3 F T F).
(4.17)

J can then be defined as the determinant of the F the deformation gradient tensor. Next, I c the first invariant of C, the right Green-Cauchy stress tensor, ψ the strain energy density function and G, the standard shear modulus. The Kirchhoff strain tensor τ s can then be written as:

τ s = GJ -2/3 F I - 1 3 I c C -1 F T . ( 4

.18)

Spatially integrate the Deviatoric Kirchhoff stress tensor Since τ s is evaluated from the fluid velocity field, it holds the fluid strain τ f as well as the solid. By integrating this updated solid strain in space one obtains the structural forces.

f = (τ s -τ f )dS. (4.19)
Here f depends on the fluid mesh as this force will be applied to the fluid of equation ( 4.1) in the last two terms.

Evaluation of inertia At this stage, it is possible to evaluate inertia as well, if needed. It requires δρ = ρ s -ρ f , the difference between solid and fluid density.

f ine = δρ(g - Du s Dt )dS. (4.20)
The inertial forces are to be added to the solid forces computed in the previous paragraph, only if δρ = 0. Dus Dt is the total derivative of the solid velocity, and g represents a possible body force per unit of undeformed volume (i.e. gravity effects).

Update the solid point coordinates Finally, the kernel positions are being updated using the solid velocity, such as:

x n+1 s = x n s + δt 2 (u n s + u n+1 s ). (4.21)
Once this step is complete, the structural forces from equation (4.19) are used to compute the error which must be below the tolerance k tolerance in order to proceed to the next time step, as defined in the flow chart 4.1.

After detailing all the steps of the FSI code used in this thesis, a test case is chosen in order to verify the accuracy against finest mesh simulation. Multiple benchmarks where achieved by [START_REF] Gil | The immersed structural potential method for haemodynamic applications[END_REF], such as free fall of a rigid and deformable cylinders, or floating or sinking balls in a closed channel, or flapping membrane under pulsating flow. There are a few tests that are commonly used to validate FSI code. For example, the movement of a free moving flap behind a cylinder has been widely studied which is a common geometry initially proposed by [START_REF] Turek | Proposal for numerical benchmarking of fluidstructure interaction between an elastic object and laminar incompressible flow[END_REF]. Yet, for this document, the benchmark used here is the deformation of a particle in a Couette flow, and it will be detailed in the next section. This numerical test has been selected because it has not been tested in the initial paper, and an analytical solution is available as well.

Benchmark case -an elastic particle in a Couette flow

Introduction

This benchmark allows us to test the validity of the physical solution provided by the numerical code. However, the numerical examples provided in [START_REF] Gil | The immersed structural potential method for haemodynamic applications[END_REF] paper are not interesting for lymphangion modelling. The flow around ellipsoidal particle in a Couette flow has been analytically solved by Jeffery [1922]. Since then, simulations has been used to extend this solution to elastic, deformable particles. When deformed by the flow, the elastic particle takes the shape of an ellipsoid. A relation between the capillary number (Ca) and the aspect ratio of the ellipse, also known as the Taylor deformation parameter (D) is shown to be linear according to [START_REF] Gao | Deformation of elastic particles in viscous shear flow[END_REF] for low capillary numbers3 . Moreover, the slope of the linear regime is shown to be one. Here, Ca 0.16 for the simulated particle, thus similar value for D is expected. Particle deformation is key element for the study of lymph as those particles can be assimilated to lymphocytes in the following sections. In addition, physical parameters of a lymphocyte can be modified by their biological functions. It is therefore important to check that the code can properly handle particle deformations.

First, the geometry and dimensions are given, then the boundary conditions are defined, finally error are estimated and compared to analytical solution.

Presentation and geometry

The geometry shown in figure 4.4 is defined in order to test the physics of the two dimensions code used in this thesis The boundary conditions are set to generate a shear flow which will then deform the particle. The boundary conditions are therefore u x = U and u y = 0 on the top wall and u x = -U and u y = 0 on the bottom wall. Initially, the particle is set as a disk, then under shear flow it will be deformed into an elliptic shape. The shape of an ellipse is given by the following equation:

1 a 2 ((x -x 0 ) cos(θ) + (y -y 0 ) sin(θ)) 2 + 1 b 2 ((x -x 0 ) sin(θ) -(y -y 0 ) cos(θ)) 2 = 1. (4.22)
Here, a define the semi-major axis of the ellipse and b defines the semi-minor axis of the ellipse as shown in figure 4.4. Furthermore, the angle θ is the angle between the long axis of the ellipse and the x-axis of the reference frame. In our case the simulation domain is set in a way that x 0 , y 0 are the coordinates of the center of the simulation box. The dimensions defined in figure 4.4 and the size of the box l = 150 µm and the initial radius of the particle r = 55 µm. In this benchmark, we set U = 1 cm s -1 . In terms of discretisation, four different fluid mesh densities are used, those are listed in table 4.1, the solid mesh density remains however the same. 4.1, the solid mesh density remains however the same.

In these four simulations a constant number of solid kernels are used and the number of fluid cells are varied. The time step (δt) is decreased along as the number of fluid computational points is increased to maintain stability of the simulation based on the CFL condition defined at the beginning of the simulation. Also, the second column of table 4.1 shows that time spend for each timestep is larger as the mesh is refined. Table 4.1: List of the discretisation and performance parameters for particle stretch in a pure shear flow. Here, time-step per second refers to the amount of time-steps achieved in one seconds. δt corresponds to the time discretisation of the simulation. nx and ny represent the number of fluid computational cells in each direction x and y. h is the length of a fluid cell, it is obtained by l nx and is mesured in centimeter. In figure 4.5, the different norm of error are shown for the x & y-velocities as well as the pressure. For example, Ω being the computational domain. i and j belong in Ω. One then can write the x-velocity L 2 error like so:

L 2 (Ω) = i,j u i,j -u fm u fm . ( 4

.23)

And for the L ∞ :

L ∞ (Ω) = max i,j u i,j -u fm u fm . ( 4

.24)

With u i,j the value of the computed velocity and u fm the finest meshs used to compare. To estimate the error for all the fields in this simulation, we compare the pressure to a zero pressure field which is a correct assumption if one assume that the particle did not start to move. Then, as the solid kernels start to move and deform the pressure solution will be different from a zero scalar field. Same reasoning applies to the y-velocity as it is set to zero at the beginning of the simulation and grow positive and negative value as the particle deforms.

Results

In order to find the proper value for the elongation of the particle one can use the variables a and b defined previously. Those long and short axis of the ellipse define the Taylor deformation parameter D, as defined in [START_REF] Jackson | A model for large deformation of an ellipsoidal droplet with interfacial tension[END_REF] but also in [START_REF] Gao | Deformation of elastic particles in viscous shear flow[END_REF].

D = a -b a + b . (4.25)
With this parameter D, we will analyse the evolution of the deformation of the particle with the capillary number Ca and compare it with the analytical solution. This dimensionless number compares the viscous effect over the elasticity.

It is defined as :

Ca = γµ fluid η solid . (4.26)
Where γ is the shear rate defined by the boundary conditions and is derived like so 2U l . µ fluid being the fluid viscosity and η solid is the shear modulus of the solid. In this case, U = 0.1 cm s -1 and the size of the simulation box is l = 15 × 10 -3 cm, so γ ≈ 133 s -1 . Furthermore, µ fluid = 1.2 × 10 -2 P for the water viscosity, and also η solid = 10 Ba. So the capillary number for these simulations is Ca ≈ 0.1596.

Also, [START_REF] Gao | Deformation of elastic particles in viscous shear flow[END_REF], defines a relation between Ca number and the deformation parameter D. In order to obtain this relation, they use the dimensionless stream function from the Stokes flow solution of a rotating cylinder. From this stream function, a velocity can be determined as well as the pressure field. At the interface between the cylinder the stress distribution is obtained from the fluid pressure field. Then the velocity, stress and pressure in the solid cylinder is then expended in the power of the Ca number. Using simulations and mathematical simplifications, it is possible to obtain the pathline of a material in function of the Ca number. Finally, equation 49 of [START_REF] Gao | Deformation of elastic particles in viscous shear flow[END_REF] writes :

D = Ca + O(Ca 2 ) (4.27)
However this equation have limits such that the capillary number must be below 0.3 and the cylinder deforation must be small. Furthermore, this paper study particle-particle interactions in a shear flow. The points of computation distribution here have been designed to take the shape of concentric circles with same ratio of surface per points. However, this way of distributing computational points is not isotropic, thus making the cloud of points with better discretisation in one direction compare to the other as shown in figure 4.6 (i). Figure 4.6 (i) show the particle discretisation at an early stage of the computation. It is slightly tilted, with θ < π 2 . Furthermore, figure 4.6 (ii) shows the interface of the particle after it reaches a plateau deformation and begins to From figure 4.7, all the mesh converge toward a similar deformation parameter D. The finest mesh represented by the simulation number 4 and its final converged value is 0.1587 approximately. Also, using equation ( 4.27) one can compare the Capillary number and the deformation parameter D. Here, the Ca number for these simulations is 0.1596 so the converged value of D should be the same. But using figure 4.7, the converged value of the deformation parameter D is lower as the mesh is finer. The finest mesh (curve 4) one can read a value of approximately 0.1587 for the deformation parameter D. This value is lower than expected but the relative error of the deformation parameter compares to the capillary number is of 0.52%.

Dispersion of particles in a lymphangion with valve

In this section, we use the 2D simulation code to the study of the transport and dispersion of a cloud of particles representing lymphocytes in a model lymphangion. Multiple parameters are varied to evaluate their influence, such as the number of particles, or the presence of valve in the flow. The flow and structuration of suspensions of deformable particles in flow is an active research field due to its complexity. Many works report on the distribution of red blood cells in a blood vessel for instance and its consequences for transport and rheology. [START_REF] Geislinger | Hydrodynamic lift of vesicles and red blood cells in flow -from fåhraeus & lindqvist to microfluidic cell sorting[END_REF], define different types of forces that apply to a particle in a Poiseuille flow and lead to transverse motion:

• Shear-gradient causes a migration of particles.

• Wall-particle interactions leads to an inward migration, toward the center line.

• Hydrodynamic interactions between particles diffuse the particles away from each others.

The first two lead to migration of particles toward the center of the channel while the last force mentioned spread the particle away from each other. [START_REF] Kaoui | Lateral migration of a two-dimensional vesicle in unbounded poiseuille flow[END_REF] investigate numerically in two-dimension the migration of a vesicle in low Reynolds number. Furthermore, non-inertial migration of vesicles in a Poiseuille flow is studied experimentally and numerically in [START_REF] Coupier | Noninertial lateral migration of vesicles in bounded poiseuille flow[END_REF]. This equilibrium state is formed maintaining the particles in a central zone in the blood vessel, as explained in [START_REF] Dmitry | Blood flow and cell-free layer in microvessels[END_REF]. [START_REF] Sylvain Losserand | Migration velocity of red blood cells in microchannels[END_REF] measured the inward migration velocity of red blood cells in microchannels and derived a universal scaling law for this migration velocity as a function of local shear rate and distance to channel walls.

When comparing blood and lymph in the collecting network, there are a few differences. Key differences are the particle nature and properties, the nature of the pulsation and the viscosity of the fluid. If one compares red blood cells and lymphocytes, they have a similar diameter. However, red blood cell have a disk shape while lymphocytes are more spherical. According to [START_REF] Hivroz | Biophysical aspects of t lymphocyte activation at the immune synapse[END_REF], a human lymphocyte's stiffness value can vary from 90 Pa to 1 MPa by order of magnitude. This wide range can be explained as lymphocytes are a family of immune cells whose activity and physiological state have a strong influence on their stiffness. On the other hand, the effective Young's modulus of RBCs (measured by Atomic Force Microscopy) varies over a narrow range between 16 kPa up to 34 kPa on average according to [START_REF] Kozlova | Morphology, membrane nanostructure and stiffness for quality assessment of packed red blood cells[END_REF]. However, it should be noted that treating cells, especially red blood cells, as homogeneous elastic particles characterized by an elastic modulus is a rough approximation compared to the standards of numerical modelling of red blood cells nowadays. Still, a simple elastic model has its own merits for a simple evaluation of cell dynamics in a flow. First, the particles are more elastic but also depend on the state of the immune system. Also, comparing the frequency of fluid pulsation between the blood system and the lymphatic system, it is lower for the lymph and the peak velocities are slower as well as defined in the introduction. Lastly, the effective viscosity of whole blood is different from lymph. However, the viscosity of blood plasma is about the same as the viscosity of lymph without lymphocytes. Both of them being body fluid, they depend on many factors but here we consider the fluid that composes lymph to be a newtonian fluid. After describing the geometry of the system and boundary conditions, we present here a set of results on the evolution of the spatial distribution of particles. The goal of such a geometry is to evaluate the effect of the presence of the valve on a pulsed flow with particle. In order to achieve this, one needs to run multiple cases with different configurations. This batch of simulations presented here is made of seven specific setups. All these simulations use the dimensions defined in table 4.2. When compared with the dimension of a real lymphangion only half a lymphangion length is simulated in order to save up computational time. It is a limit of this study to not have proper sized lymphangion, as it would contribute to dispersion of particle to a greater extent. Particles are chosen to fit the lymphocyte size reported in literature. As the number of particles varies with the immune response, it is interesting to test different particle concentrations. In this study, we focus on three cases with respectively 1, 36 or 72 particles, which corresponds to a density of approximately 2, 80 and 160 particles/mm 2 respectively. Furthermore, it is important to notice that the initial position of particles in the simulation box is random. For this geometry, a set of seven simulations is summarized in tables 4.3 and 4.2.

Presentation of the geometry

In terms of dimensionless numbers the Reynolds, the Womersley and the Strouhal numbers are defined as:

Re = vL ν , Wo = L ω ν and, St = ωL v (4.28)
With characteristic length L = 0.3 mm the diameter of the vessel, the kinematic viscosity ν = 0.8 mm 2 s -1 , the characteristic velocity v = 1.5 mm s -1 and ω = 0.33 Hz for the pulsation frequency. The Reynolds number correspond to the ratio between inertial forces and viscous forces, while the Womersley number compares the ratio of inertial pulsatile flow effects to viscous effects. Here, the Reynolds number is approximately 0.56 and the Womersley number is 0.18 approximately. Therefore, the unsteadiness of the pulsatile flow is not to be considered as the Wo number is below one. As the Reynolds number is small, it is reasonable to assume that the inlet boundary condition , far from valves, can be considered as a Poiseuille velocity profile whose amplitude instantaneously follows the imposed pulsed flow described below The Strouhal number is 0.066 and it compares the ratio of velocity between oscillatory velocity to ambient flow velocity.

Boundary conditions

In order to simulate what would happen in a long chain of lymphangions, we apply periodic boundary conditions on the velocity and particles. In other words, the velocity profile of the right border is used as the velocity profile on the left border. With regard to the particles, it allows them to leave from the right side of the simulation box to re-enter from the left side. A multiplicative coefficient is used to modulate the velocity, and reproduce a pulsed flow. From the immersed boundary condition method, the particle position modification was simple to implement. In fact, a particle position is defined by the average position of all its solid kernel position. Then to move a particle from the right boundary to the left, one can simply subtract the length of the simulation box to the x-coordinate of each solid kernel position. So its orientation relative to the simulation box and its size are conserved, since this operation is a geometrical translation. Figure 4.9 shows the evolution of the boundary condition applied to the left wall. This periodic boundary condition cycles with a period of 2 s. For a short period of time, the flow is reversed, as imposed by the boundary conditions. The flow reversal due to this pulsatile flow needs to be implemented in order to take into account its effect of valves, which in turn periodically modifies the geometry and therefore influences particle dynamics.

t(s)

u in (cm s -1 )

t n t n + 0.5 t n + 2 -0.05 0 .1
Figure 4.9: Periodic boundary conditions in time for the inlet x-velocity imposed on the left boundary condition. From a mathematical point of view the function periodicity can be written: uin(t) = uin(t + 2). Reader should note that uin is the maximum velocity of the Poiseuille profile.

As defined by the Poiseuille flow, v in , the velocity in the y-direction is set to zero. The top and bottom walls are considered as non-slip boundary conditions. With the velocity equal to zero at all times.

Mesh convergence study

In order to estimate the error produced by the code, a mesh convergence study is achieved. Since there are two materials, namely the fluid and the solid to discretise, only the fluid mesh size is varied in this study.For the solid structures, no meshes are available, only a collection of solid kernel points. However, the density of solid kernel points used to define the solid does modify its discretisation. For this convergence study, only the valve is considered, particles are removed of the simulation. Leaflets are made of 24300 solid kernel points each, so a valve is supported by 48600 kernel points. In order to produce a leaflet, a rectangle is first produce of 270 times 90 solid kernels, then this rectangle is deformed via geometric transformation into an elongated shaped leaflet. Nevertheless, the simulation time increases with the amount of solid kernel used. So apart from this mesh study, for later simulations a balance has been found between computational accuracy and time. Fluid mesh density refers to the amount of pressure computational points divided by the surface covered. Similarly, solid kernel density, is defined as the amount of kernels divided by the initial surface covered by the solid. A condition linking the fluid mesh density and the solid kernel density is provided in [START_REF] Gil | The immersed structural potential method for haemodynamic applications[END_REF], which should be at least 2 solid kernels for a single fluid computational point, in each direction. In fact, they tried different types of kernels but only the one defined in [START_REF] Charles S Peskin | Flow patterns around heart valves: A numerical method[END_REF], equation (6.27), is properly used and detailed. So for issues of clarity and stability, the Peskin kernel is the one used in this section. Table 4.4 shows the list of all meshes chosen for the mesh convergence study. Next, relative errors for the velocity fields and pressures are plotted in figure 4.10. In figure 4.10, the L 2 and L ∞ -norms are displayed in function of the x mesh size, called h. As no analytical solution exists, we compare the whole velocity and pressure fields solution to the finest mesh. Table 4.4: List of the discretisation and performance parameters for the valve study in a straight channel. Here time-step per second refers to the amount of time-steps achieved in one seconds. δt correspond to the time discretisation of the simulation. nx and ny represent the number of fluid computational cells in each direction x and y. h is the length of a fluid cell, it is obtained by l nx and is mesured in centimeter. The error for the velocity in the x-direction, y-direction as well as the pressure field are compared to the finest mesh case. This case is the run number 4 from table 4.1. Orange and yellow linear slopes are here to represent respectively the error in O(h) and O(h 2 ).

L ∞ norm has a higher rate of convergence than what can be observed in [START_REF] Gil | The immersed structural potential method for haemodynamic applications[END_REF] as it is not impacted by immersed boundary additive term called f in equation (4.1). As a matter of fact, the pressure field is only estimated locally at the computational point x = 0.03 cm and y = 0.015 cm. The reader can appreciate the proportions as well as the size of particles and how they spread in the channel. One can observe that the particle tends to gather in the center and distribute themselves properly across all the entire channel. Also, for capture time 22.5 s, the Poiseuille flow is going on the opposite direction, since the input velocity frequency is 0.5 Hz as shown in figure 4.9 as u in is a periodic function. This frequency is higher compare the contraction frequency showed in the lumped model chapter, because the simulation does not represent a complete lymphangion. Here, in figure 4.12 display the entire simulation domain with velocity fields in black arrows and the valve and particles in red. At time 2.5 s the flow is reversed, and it is possible to see how the valve is impacted by the flow reversal. As these simulations with valve are shorter than ones without valve, the spreading of particle in the y-direction is not complete. Figures 4.13 & 4.14 characterize the migration of particle towards the center distribution of particles at specific times Particles are spread randomly in the simulation at time t = 0 s. These figures also show that there are no large clusters of particles forming and that they tend to spread along the x-axis.

Study of the distribution of particles

(i) (ii) (iii)
(i) (ii) (iii)
In figure 4.13, the evolution of the distribution of particles along the y-axis is shown, with a time step of 7.5 s and by counting particles in layers of 0.3 mm width in the y-direction. Comparing between different particles densities in the valveless cases, figure 4.13 shows that for all simulations including simulation C, the concentration tends to amass in the center. At t = 7.5 s, there are approximately six particles per range between 0.6 and 2.4 mm. And at t = 30 s, the central range has 15 particles, less than half the total amount is concentrated in one tenth of the channel width. It is not the case for simulation D, but a similar process occurs: the number of particles decreases near the walls, increases in the center. But it is not possible to observe a single central channel, the three main ranges have 16 particles each, although the distribution is not as narrow as in the previous case. In this simulation, more than half of the particles are spread in less than a third of the channel width. This is likely due to an enhanced shear-induced diffusion due to a greater particle concentration and also possibly to the fact that the configuration has not reach steady state.

In the case of valve simulation (case F & G), the maximum time is only 5 s, so stability is not reached yet as it can be seen from figure 4.18. So the results of figure 4.14 are not converged into a stable state. However, similar to the cases without valve studied in the paragraph above, the particle concentration is increasing in the center of the channel and decreasing in the region near the wall. In terms of the distribution at the latest time step, it is similar to the simulation C case, i.e. no homogeneous front appears. The part linked to the presence of an obstacle in the flow or simulations has not reached a time independent state and thus this state cannot be evaluated. tions without valve. The pulsatile nature of flow does not have a large impact on the average y-position over time for valveless simulations as no trace of the 2 s frequency can be observed. While the average y-position of simulations B and C tends to quickly reach the center of the channel, defined at 0.15 mm, the D simulation shows a slightly different behaviour. This simulation has 72 particles, and thus has the highest density of all valveless simulations. The initial positions of the particles are generated randomly. As a consequence, in the D simulation, a cluster has formed and will take more time to diffuse properly than the total simulation time. Nonetheless, the final average y-position of all particles is close to the center of the channel in all valveless cases. Figure 4.15 graph (ii) represents the average position of all particles on the y-axis, for simulations with valves. The total simulation time with valves is 5 s while it is 30 s in simulations without valves. For simulation E with only one particle, effects of the valve are quite important and the cyclic variations are large. For the F simulation with 36 particles, the spreading action on the x-axis is more consequent and the cyclic variations are smaller. Lastly with the G simulation, there are a small variations and the overall value is close to the center of the channel.

It is possible to compute migration velocities based on the initial slope for figure 4.15 (i) and the average slope for figure (ii) where simulations have a duration of 5 s. Note however that, simulations D and G start with a cloud of particles that is already more or less centered, and show mainly fluctuations around this position. For simulation B without valve and only one particle one obtains a migration velocity of -0.61 µm s -1 , using the initial slope from 5 s to 15 s. Furthermore, for the simulation C the velocity almost doubles from the B case and reaches -1.18 µm s -1 . In the valve case, the situation is the opposite. E simulation has a velocity of -1.63 µm s -1 and the F one has velocity of 0.65 µm s -1 .

The difference between the two velocities for a single particle case can be attributed to the presence of the valve. The velocity of migration depends on shear that the particle is submitted to, so the initial position does impact the velocity as well as the presence of valve. However, the difference between the two velocities for 36 particles is less obvious as it might come from an equilibrium between the diffusive force spreading particles in a repulsive way from the center to the walls of the channel and migration force pushing particle toward the center. 4.16 shows the distance in the y-direction between the walls and its closest particle. This is a measure of the thickness of the particle-free layer near walls. (i) represents only valveless simulations, and (ii) simulations with valve. In the valveless cases, the observation formulated in the previous paragraph are confirmed. After 8 seconds of simulation D, values plateau or oscillate around same value. Furthermore, the distance from the top and bottom wall are similar. For simulation C, such phenomenon occurs as well but after a long period. Both curves, C t & C b raise for 18 seconds, here subscript t means top and b means bottom. However, the distance from top wall falls back while the bottom one keeps rising. In figure 4.16 (ii) the presence of the valve causes large peaks. The average migration velocity is faster when there are valves in the lymphangion. The overall displacement is bigger after 5 s in 4.16 (ii) than after 5 s in 4.16 (i). Due to the reduction in channel width, shear rate is higher in the valve section, which enhances the migration towards the center. However since the particles are at different distance from wall, due to the Poiseuille nature of the velocity profile, it means they also travel at different speeds. In fact, a particle closer to the wall is advected slowly compare to a particle travelling in the center of the channel. So, they do not get close to the valve at the same time nor at the same frequency.

F b Ft G b Gt (i) (ii)
The distribution of particles as well as the velocity field with arrows are represented in figures 4.17 & 4.18 at specific time. For the valveless simulations, capture at 22.5 s shows a reverse flow configuration, while the simulations with valve happen for the capture at 2.5 s. It is possible to observe from figure 4.17 that with fewer particles the y-displacements are limited. In fact, the graph (i) shows a very few changes in direction and most of the curves are straight. While in the graph (ii), curves tend to oscillate more and overall travel larger distances in the y-direction.

Multiple ideas can justify such phenomenon such that the diffusive force that spread the particles can be affected by the density of particle. Furthermore, it could be linked to the initial position of particles, but no statistical simulations study were run. The particle density however does not seem to influence the movement of particles in figure 4.18 for simulation with valve. As a matter of fact, most of the displacement of particles in the y-direction are linked to the valve. It pushes particles toward the center then draws them back toward the walls after. Such movement produces larger amplitude in the particle displacement on the y-direction. In order to obtain the different migration velocities several y-ranges are defined, similarly to distribution of particle study. So all the particles in a band is used to compute the average in space. The time average is done in such way that a particle can only be counted once, but across different y-ranges. The migration velocities are low compare to the advection velocity, by one order of magnitude roughly. time progress. This migration toward the center is no depending on the upper or lower side of the lymphangion. Also, the migration velocity reduces with time in both cases with 36 or 72 particles. However, the particle density appears to affect the intensity of the migration velocity. It can be linked to the effective viscosity modified by the increased particle density. Such hypothesis should be confirmed by the cases with valve. Furthermore, the symmetry between the top and bottom region remains the same. In figure 4.13, the two cases have valve. It has been observed in the previous subsection that the valve contributes to the mixing of particles with figure 4.18 for example. In the simulation F, the average velocity for the particle below the center line have negative velocities. This behaviour can be linked to a group of particles moving toward the wall after passing around it. The simulation G does not show this behavior because the density is larger and the averaging process smooth these local phenomena.

Velocity profiles

In this subsection, a closer look will be taken at the velocity profiles of the lymph in a lymphangion. The Reynolds and the Wormersley number are low as computed in Presentation of the geometry. Therefore, the fluid inertia is low compared to viscous effects as well as the influence of the pulsatile flow. However, secondary lymphangions have valves which perturbed the flow. Consequently, it is possible to observe such impact on the width of the established Poiseuille solution. First, simulations C & D will be considered in figure 4.21, which are valveless. Then the simulation F & G that present a valve will be observed in figure 4.22. smooth whereas the other profiles are rougher. This suggests that the flow achieves a steady state after at least 21 s and is only perturbed by the particles. Their position is defined randomly at time t = 0 s, and most of the perturbation are present near y-position 0.1 cm and 0.2 cm. These two bands are equally distant from the center, and does depend slightly on the particle density. Furthermore, the intensity of these perturbations is larger in the simulation D.

For simulations with valve, in figure 4.22, profiles plotted for x values of 0.06 cm, 0.09 cm & 0.12 cm the impact of the valve is observable. In fact, negative values can be seen near the walls, because a recirculation zone develops between the wall and the valve. However, after 3 s, the valve lie over the wall, as can be seen on the plot x = 0.12 cm. Before 3 s, the valve is not open and the parabolic profile is confined between y = 0.006 cm & y = 0.024 cm which concentrate the particles and produce perturbations on the parabolic profile.

So particles do not modify to a great extent the parabolic profile, however the perturbation are concentrated at the same position. The valve however impact the flow and tend raise locally the particle density which generates more perturbation.

Position of the valve

For this subsection, the valve displacement is the main focus. The initial setup is described in figure 4.8. The boundary condition are identical to the previous subsection, a Poiseuille flow for the inlet and outlet, and non-slip wall for the top and bottom simulation boundaries. For the valve, they are embedded to the wall, and the root of the leaflets do not move during the whole simulation. The two-dimensional simulation used here contains the two leaflets of a lymphangion, but only the tip of the lower one will be discussed. Simulations are run with 5 seconds of physical time, which provides more than 2 contraction cycles. In figure 4.23, the initial position of the leaflet tip is defined at the (I) point, then opens up to reach the point (II). At this moment the flow has changed its direction and instead of closing the valves it keeps on opening but the tips retracts (II → III). Then the flow returns to a forward direction and keeps opening the valve also, the number of particles influence the valve tip as it can be observed in figure 4.24. (III → IV). At (IV) the valve has reached its full extension and cannot extend more. It will retract to (V) position where another reverse flow event occurs, and then into (VI) where a similar event also happens. It is interesting to notice that the valve keeps on opening even with reverse flow. This might be linked to the reverse flow not being strong enough, or the reverse phase having a period too short. The initial opening phase pushed the valve from an opening radius of 20 µm up to approximately 78 µm. Then the second phase reached an opening radius of 100 µm. However, the third phase does not increase the opening radius by a significant margin. A normal lymphatic valve have an approximate value comparable to the one computed here, and should return to a closed state when backflow occurs.

With regard to the different amount of particles in the simulation, it alters the opening movement of the valve, but it can be considered negligible. However, the path taken by the valve varies slightly between simulations. The largest modification comes from the G simulation which has the largest amount of particles. It is also interesting to notice that whatever the path takes, the loops produced by the reverse flow occur at the same y-coordinates. So the particles in this case mostly affect the elongation of the valve rather than its opening radius. This effect can be attributed to multiple causes, but the main one would be that the valve is more free in the two-dimensional case than in real life. In fact, [START_REF] Li | The effects of valve leaflet mechanics on lymphatic pumping assessed using numerical simulations[END_REF] used a variable thickness (Young Modulus) of alongside the valve main axis to approach a more natural behaviour. Such setup has been tried as well during this thesis, but is not presented in this manuscript. After showing results of simulation with pulsed flow, a comparison with regular Poiseuille flow is considered. The initial position of particles is different from simulations A to G, but after a period of time it does not influence the results.

After comparing different particle densities and valve impact on the lymphatic flow, a next step is to observe the influence of wall. Since those actively contribute to the pumping aspect of the lymph displacement, numerical simulations will be performed to understand this phenomenon. However, it is still subject to debate, and no geometry properly captures the whole physics. The following section aims at remaining simple and only focus on 2D geometry.

Towards simulation of pumping dynamics

Similarly, a lymphangion contracts its wall to move lymph, we study the effects of the wall on the displacement of the valves. Furthermore, the flow variation in the channel modifies the position of the valve relative to the wall allowing the fluid to pass or else be stopped. We compare the results in input/output flow rate, and different geometries. simulation domain total height 0.5 mm total length (L) 1.5 mm movable wall thickness 0.02 mm distance between walls (l) 0.3 mm leaflet width 0.03 mm length ≈0.33 mm Table 4.5: List of the main physical values for the wall impact on lymph flow in a lymphangion test. In table 4.5, the values of the geometry used for this section 's simulations are summed up. The length is defined based on the half lymphangion geometry from the previous simulations. The height is defined to allow the movement of the wall on the y-direction. Leaflets and walls have been chosen based on the geometry presented by [START_REF] Daniel | Integrated geometric and mechanical analysis of an image-based lymphatic valve[END_REF] and Zawieja [2009]. The geometry proposed here imitates the experimental setup presented in appendix. Its main objective is to move the walls as defined in figure 4.25 to compress or expand the lymphatic volume. In order to achieve this goal, multiple strategies were used.

Presentation of the geometry

Methods

Three approaches of this study were used to simulate effects of contraction on the wall's displacement. Boundary conditions are set such that:

• Flow rate boundary function is forced on fluid field,

• External displacement function is imposed on the solid walls,

• Additive force function is enjoined on the solid walls.

Initially, it was planned to use boundary conditions to impose a flow rate on the top and bottom wall to force the fluid to produce a wall displacement. However, a boundary condition function along the y-axis must be determined to adjust. Moreover, the link between the flow rate and the pressure gradient exerted by the fluid on the wall can be determined but the one between the pressure gradient and the wall displacement seems to be non-linear. Therefore, this approach is easier to implement numerically but harder to adjust physically.

Then it was decided to force the displacement of wall. This setup was not tried since a strong control over the wall displacement is impossible in this configuration. It can be achieved by adding the wanted displacement in the equation (4.21). However, since this step is computed after the forces, the interaction between fluid and solid is not properly achieved. Also since the physics is altered from the solid point of view, the convergence study did not succeed. And it was then not kept for any of the following results.

Lastly, a force field was applied to the wall to simulate the contraction and relaxation of the lymphangion. Compression of lymph worked properly, but the inflating of the lymphangion gave rise to instabilities. Those instabilities are known in fluid mechanics as Rayleigh-Taylor instability. So the wall was behaving as the interface between fluid of different density even if in our case the two fluids had same density, but with a force applied to the interface. This method is easier to implement and also provide more interesting results. Here the wall are assimilated to an interface, and the contraction motion is achieved in a way that resemble the biological way as if muscle cells were stimulated.

The following part focuses on the mesh study, those are run for all different geometry cases presented in this section. Doing so allows the mesh to be verified, as well as the physics. In fact, the pressure and velocities fields are verified against the finest mesh simulation.

Mesh convergence study

Similarly to the previous section with particles flowing in Poiseuille flow, a mesh convergence is achieved. The list of simulations runs for this study is shown in table 4.6. Furthermore, the figure 4.26 shows norm L 1 , L 2 and L ∞ , long side h which represents the mesh size in the x-direction. In both cases, the parameters or the solution are tuned to the finest simulation computed. The first case appears to be diverging, the norms of the error are above one for both norm L 1 , L 2 . However, the three other cases present a standard slope. Similarly to previous setup (without movable walls), the pressure has low error for the latest case. Furthermore, it appears to maintain a constant slope while the velocity does not. These results were obtained using the latest setup simulations and with a physical time of 0.8 s.

Details on different geometries to model pumping

The main concept is to observe how the wall compressing or expanding the fluid could impact its flow. Multiple approaches to this problem have been considered, and only one was properly tried. The initial one was to apply Dirichlet type of boundary condition to enforce wall movement. However, a more direct approach was retained which consists of applying the force field directly onto the wall. This force takes the form of an additive constant to the force field from equation (4.19). From this setup two different actions of a lymphangion can be defined:

• Push configuration -contraction motion -where the fluid is compressed in between the walls and is propelled outwards through the left and right boundaries.

• Pull configuration -relaxation motion -here, the fluid is absorbed by the lymphangion, so it moves inwards through the left and right boundary of the simulation box.

For the initial test, it has been decided to have no solid attachment points of the wall relative to the simulation domain. In terms of boundary conditions a

Neumann of zero are applied to all the walls orthogonally. In order to maintain the mass balance, the right wall is still defined as the inlet or outlet, depending on the configuration. This is achieved by a multiplicative constant applied to the boundary condition of the right wall. Overall, all this work on the boundary condition helps comply with any trends the wall can take. For example in figure 4.27, one can observe that the whole domain is pushing fluid to the left direction. This shows the requirement for a boundary condition to hold the wall in position and not be advected by the fluid. Furthermore, it is possible to see the effect of a pull configuration in figure 4.28. Here the fluid enters from the left and right walls of the simulation box. The valve is still in the initial position state, as the velocity field has a small magnitude in the center of the simulation. In this type of configuration, the movable walls are advected up and down meaning the lymphangion radius increases. However, the movement should be caused by the fluid forcing the inflation of the container and not the wall itself. This was not tries as it needs to redesign the whole simulation and treat the boundary conditions properly. Lastly, figure 4.29 represents the second geometry. It was only tried in a push condition because of the reasons mentioned above. Furthermore, the Young's modulus of the wall and valve are modified between the two geometries. Here to fix the advection of the wall, a small square near the intersection between the valve and the wall forces the x-velocity to zero. These patches only allow the walls to move up and down, setting the solid part in the whole simulation box. As this is a push configuration, fluid should flow out of the simulation to the right and left sides. Compared to figure 4.27, the results are more physical. Furthermore, outside the lymphangion on the left boundary, the fluid enters the simulation domain to press the lymphangion as the walls have additional force. On the right boundary the magnitude of the fluid is large inside the lymphangion than outside.

Observing such behaviour allows us to quantify the flow leaving the simulation from the right and the left in order to estimate difference of flow provided by this geometry.

Flow analysis

Here we are comparing the evolution of the flow in a push configuration. The flow is estimated at the left and right walls of the simulation between the top and the bottom walls, as the flow of the lymphangion. It is important to notice that in this case both before and after the valve are contracting in the same time. As showed in the lumped model chapter, that the contraction only accounts for 10% of the whole period cycle. Furthermore, two lymphangions contracting at the same time can occur in a network it has been mentioned that they tend to synchronize. In such way, if one is contracting the other one happen to be in a relaxing phase. Figures 4.30 & 4.31 show the flows and diameter after and before the valve. Since the simulations take a long time to complete, the time frame studied here is only the initial phase of the contraction. However, it is possible to observe how the forward motion propagates. The initial setup should lead to equal propagation of fluid to the right and the left, but the valve forces the left-hand flow to reduce in intensity. On the graph 4.30, after approximately 0.15 s, the right flow keeps on growing while the left flow keeps on decreasing. A sum of the two represents the overall valve impact. However, these results should be considered with care as they only represent the initial phase of the contraction. The diameter for example on figure 4.31, should keep decreasing, after 0.11 s for the left wall and 0.18 s the walls expend again. This phenomenon can be linked to an equilibrium between the diameter and maximise the flow.

Chapter 5

Conclusion

Summary of the thesis achievements

This thesis provides a study on a multiscale approach of the secondary lymphatic network. The first scale considered a chain of lymphangion forming a vessel using lumped models. It provided insight on the pumping function of a lymphangion as well as its synchronizing ability in a chain to produce a forward flow in a vessel. This might lead to interesting findings in order to improve condition of people with light lymphoedema symptoms. As of today there are little to none clinical treatments, and the most effective one is compression garment. Lumped model is the standard for computational lymph network but refining models and computational methods could lead to better understanding of lymphatic flows, including in its pathological aspects. As there is no refined lymphatic network simulations in the literature, one could find or proper medication except compression garment. The second scale studies part of a lymphangion with coupled fluid and solid interaction in two-dimensional simulation. With these FSI simulations, we compared various cases to find whether presence of valve in the lymphatic flow or the density of particle would have an impact on the transport capacity. Furthermore, different geometries have been tried, with or without movable walls, which allowed us to confirm that secondary lymphangions only actively squeeze fluid and passively inflate.

In order to produce these results, we studied a lumped model approach, and extended it to a fully integrated approach, using the latest findings. In terms of modeling, a recent constitutive equation was used, including the complete set of ODEs for the valve behaviour as well as the contractile effect. Using results of in vitro experimental studies on rat lymphangion vessels, from Davis et al. [2012] and [START_REF] Bertram | Consequences of intravascular lymphatic valve properties: a study 135 of contraction timing in a multi-lymphangion model[END_REF], the model parameters were adjusted. This allowed to understand the influence of parameters on the different phases of contraction, relaxation of the lymphangion cycle, as well as the relation between the pressure and diameter. From this initial study, a simple confluence and bifurcation of lymphatic channels were investigated to observe how different vessels would react to different pressure loads. As shown by [START_REF] Philibert | Anatomie, physiologie, pathologie des vaisseaux lymphatiques considérés chez l'homme et les vertébrés[END_REF], lymphatic vessels have different lengths and various number of lymphangions per vessel. This observation led to a study to check if there is an optimum value of lymphangion 105 number per vessel. The results were observed in terms of work per lymphangion and frequency of contraction. Then we studied a larger, complex network by digitizing the geometry of a real leg skin secondary lymphatic network. This secondary network posseses more confluences and bifurcations than what has been studied in the literature. It led to interesting results on the synchronicity of the lymphangions depending on the pressure difference between the inlet and outlet of a vessel and thus its position relative to the network. Three categories were made to classify these different behaviours :

• chaotic: here no synchronisation were found between the lymphangion in the vessel,

• mixed: Intermittent synchronisation and a wave-like pattern was observed,

• synchronized: the lymphangions were all in-sync and produced wave repeatedly.

This classification might give rise to threshold value based on pressure differences between the inlet and outlet or even external pressure. As a matter of fact, compression garments are already in use, but knowing the secondary network structure of a specific person key area might require higher or lower external pressure to balance lymph flow and thus avoid fibrosis. For the twodimensional approach, two geometries are used to model a lymphangion. The first uses a straight channel with valves and particles, the second one includes the movable walls. For the first geometry, different observations were made on the valve movement, the particle distribution and on the influence given by a pulsed flow. In fact, in fluid mechanics, a Poiseuille flow is considered timeindependent but this is not completely true for lymph flow, and thus would have an impact on the particle distribution. Results on the valve movement were mitigated since the valve did not manage to close itself under the reverse flow phase. However, the particle distribution is modified by the valve presence, and they tend to gather in the center quicker than without the valve. The distribution of particles in a straight channel and in a steady flow has been studied in the case of red blood cells, such as presented in [START_REF] Sylvain Losserand | Migration velocity of red blood cells in microchannels[END_REF]. Two phenomenona enter in competition, one where the particles diffuse from the initial configuration. This spread happens as the particles are advected, so in our case in the y-direction. The second is a force that attracts particles toward the center of the channel often mentioned as "lift". These forces are observed in the simulations performed here. Also, the presence of the valve modifies the equilibrium between the two phenomena. Regarding the valve simulations, the stiffness and the geometry were not properly designed, as they did not close properly. The movable wall simulation allowed us to observe the flow difference upstream and downstream of the valve. Two different configurations emerged from this geometry, one where the fluid is pushed away from the valve, and another one where the fluid is pulled toward the valve. From the first configuration we were able to observe the flow before and after the valve to find that the valve shape allowed the fluid to move in the forward direction. The second configuration did not lead to interesting results.

Discussion on the lumped model chapter

When we compare diverse experimental results to the results from the simulation, the number of parameters are excessively large and concurrent on same physical phenomena. One good example is the ratio of time spent between the contraction and relaxation time, it has three parameters of adjustment.

The same occurs for the amplitude of contraction, with five other parameters.

Even if this problem allows multiple solutions, a simplification of the constitutive equation is welcome. Also, the ODE used for the contraction/relaxation system could benefit from a simplification of the constants used, even if the frequency response might be difficult to then integrate to the model. The use of artificial intelligence could be interesting, in order to fit the large number of parameters. The simple network section with unique bifurcation or confluence has delivered fascinating results, mainly on the frequency of contraction as well as on the work per lymphangion. In the case of a confluence, the afferent vessel is more active than the efferent ones and has on average the same contraction frequency. For a bifurcation, we demonstrate quite the opposite with the activity being similar between afferent and efferent vessels, but the contraction being on average larger for afferent vessels. On the network simulations, we used two different drawings of the lymphatic leg skin network. In terms of the overall shape, the first network proposition has reasonable amount of inlet and outlet properly spread in the leg tissue. However, for the second network, the number of inlets situated on the top part is quite significant and led to vessels with little bifurcation and confluences, and little interaction between each other. Furthermore, the simulations used eighty seconds, and it would be entertaining to allow them run for a longer period of time. This would allow the study of average inlet and outlet flow as well as compare different sides of the network to determine whether major routes are preferred or if the load is distributed across all vessels.

Discussion on 2D simulations chapter

To begin with, we described the immersed structural potential method employed to solve the numerical geometries used in this section. This method worked, but some improvement could be made on the efficiency of the solid part. As a matter of fact, the fluid section are vectorised and in addition, it used an FFT to solve the Poisson equation making the fluid solver fast. However, the solid section was not vectorised, and MATLAB spent time doing the interpolation and solving the solid. After verifying the validity against analytical solutions, different geometries were used. Initially, pulsed flow in a straight channel is used to model the dynamic of a lymphangion. Multiple topics were studied such as the spread of particles with or without valve, and also depending on the density. Also, the movement of the valve was explored, but no conclusion could be achieved from these simulations. Variation of the Young's modulus has been considered, as well as modifying the thickness gradually. These modifications were tried in order to stiffen the base of the valve, in order to imitate the three-dimensional valve behaviours.

At that point in time, another geometry was tried which took the wall motion into account however the pull configuration did not provide any results. As a matter of fact, the method of exerting forces to the wall is not how lymphangion works, such method could not work, and only led to instability in the simulation. A lymphangion is filled because the fluid is squeezed in from the previous lymphangion and not aspirated in by the lymphangion itself. So in order to get the pull configuration to work one needs to impose a Dirichlet velocity at the lymphangion inlet. Then, fill it until the walls are extended enough (estimated by the constitutive equation from the lumped model). Only under those circumstances it would be possible to shift paradigm and enter a push configuration where the wall actually contracts and a force can be applied to the wall.

Future work

With regard to the lumped model chapter, there are two main axes of research that could bring interesting results:

• Improve the constitutive model and include the chemical components.

• Tune lymphangions parameters for the network simulations. This constitutive model needs to be worked backward. In other words, it needs to be built from experiments results instead of analytically. However, such procedure requires lots of rigorous experimental work as most of lymphangions come from rodent and we are trying to dimension the simulations to a human lymphatic network. For example in [START_REF] Windberger | Whole blood viscosity, plasma viscosity and erythrocyte aggregation in nine mammalian species: Reference values and comparison of data[END_REF], sample of blood are used from different species of mammals and human seems to stand out. Even if, using rodent data can help to simulate a whole human lymphatic network one will require more data from human samples. For the network simulations, it would be interesting to play with the external pressure to observe the effect of compression garment. Also observe whether an increased viscosity would reduce or increase the pumping efficiency of a simple network. Furthermore, from Sappey drawing in figure 1.2, focus on bifurcations and confluences and compare to the results obtained in our simulations of three-vessel networks.

The direct simulation of fluid-structure interactions in a lymphangion that was performed using 2D numerical modeling in this work could benefit from improvement in several directions:

• Update toward three-dimension simulations.

• Include chemical response in the wall from the concentration model.

• Improve the lymphocyte model.

Going toward a three-dimension code would help to model the valve mechanism properly. The chemical response has already been included in two dimension by [START_REF] Li | The effects of valve leaflet mechanics on lymphatic pumping assessed using numerical simulations[END_REF] and shown decent results. The lymphocyte model does not seem to have been included in the numerical simulation of lymphangion. Furthermore, in term of geometry maybe use a three lymphangions vessel to not be disturbed by the boundary conditions effect.

Conclusion

The approximation given by the static beam theory gives coherent values. However, the maximum displacement is a bit too small compared to what one should be seeking for the set up presented here. One way to get closer would be to reduce the physical property of PDMS to get E 10 times smaller, or get thinner walls or valves to allow larger deformation.

Pressure -Pressure drop in low Reynolds

Here, the Reynolds number is low so the slope of the pressure in this case can be approximated to 16 Re . The objective of this part is to demonstrate the Darcy charge factor for this set up: a laminar flow in a two-dimension tube. As shown in figure B.1, h is the height between the two walls and L is the considered length of the flow.

As we proved in the previous part, the pressure only depends on the x-direction variable such that:

dp dx = 2µu x (y) y(h -y) . (B.12)
Then, the pressure difference between the inlet and the outlet can be obtained by integrating the previous equation:

∆p = 2µu x (y)L y(h -y) . (B.13)
Furthermore, the charge factor can be written using the pressure difference and other parameters written by Darcy and improved by Weisbach such that:

f D = ∆p2D h ρLU 2 . (B.14)
With D h the hydraulic diameter, L the length of the flow in the x-axis and U the fluid velocity. By replacing (B.13) in equation (B.14), one obtains the following formula:

f D = 4µu x (y)LD h y(h -y)ρLU 2 . (B.15)
In our case, the hydraulic diameter is four times the cross-sectional area of the flow divided by the wetted perimeter of the cross-section, in mathematical terms:

D h = 4 h 2 4h
. Then by replacing y by h/2 (and u x (h/2) = U ), one can write:

f D = 4µU Lh h 2 (h - h 2 )ρLU 2 . (B.16)
Then by simplifying the previous equation, we get: By integrating the velocity over the x-axis, considering with the y-velocity being 0 at the wall, so when x = r and x = -r, one can obtain: So from this relation, it is possible to deduce the parabolic profile of the pressure along the y-axis of a lymphangion in the case of a movable "rigid" walls.

f D = 16µ hρU . (B.17
V a = -
Their position was extracted from Sappey [1874] drawing by a software developed by [START_REF] Mitchell | markummitchell/engauge-digitizer: Nonrelease[END_REF]. A vessel number is then defined to identify all vessels but they are not used in the data set. This use the fact that intersection of the secondary network only have a maximum of three vessels. A vessel number can be characterised by the link between two connections. Table C

List of Figures

1.1 The graphs from top to bottom represent the contraction memory (i), the contraction intensity (ii), the pressure from contraction (iii) and the overall pressure (iv). Here, the contraction memory change from 0 (contracting) to 1 (relaxing). This allows the contraction memory to switch from contraction phase to relaxation phase, which then influences the contraction pressure. List of Tables

  Figure 1.1: (i) Picture representing lymph and blood circulatory system and how they connect with each other. (ii) Inner leg lymphatic skin network from large drawing of Sappey [1874]. (iii) Microscopic image showing the profile section of a lymphangion with the two leaflets from Moore and Bertram [2018].

Figure 1 . 2 :

 12 Figure 1.2: Lymphatic vessels of the leg skin (Vaisseaux lymphatiques superficiels ou cutanés du membre inférieur in french) this illustration can be found in Sappey [1874].

  Figure 1.3: Schematic representation of one lymphangions with its neighboring ones. The Diamter here represent the lymphangion diameter, and the length shows the overall length of a lymphangion. The Valves are shown in two different state. They are oriented in a direction allowing it to close if the pressure difference moves the fluid backward, and close the valve otherwise if it moves forward.

Figure 1 . 4 :

 14 Figure 1.4: Summary of elasticity values (Young's modulus) for different immune cells, from Bufi et al. [2015].

Figure 2 . 2 :

 22 Figure 2.2: Different types of valve studied from a numerical bio engineering point of view. (i) Aortic valve for numerical simulation from Spühler et al. [2018], (ii) idem from Loon, van [2005], (iii) Lymphatic valve from experimental result from Watson et al. [2017].

Figure 3 . 1 :

 31 Figure 3.1: Sketch of a lymphangion's side view, exposing the valve and the container part. (i) The position of different pressure points and diameters are provided. As well as the different flux between the pressure points (ii) The short arrows represent the wall motion of a contracting ith lymphangion. The short arrows in the i + 1th lymphangion show the valve opening. Large arrows display the different pressures a lymphangion is submitted to.

Figure 3 .

 3 Figure 3.1 represents a lymphangion with computation points for pressure and diameter. Three pressure points are estimated directly: • P v,[i] the pressure right after the valve,

  equation, k open and k close are the opening and closing compliances.

Figure 3 . 3 :

 33 Figure 3.3: Evolution of the valve state in function of time. The left graph is the closing process. The right graph is the opening process. The time scales are not relevant as they are dependent on δP , and the opening or closing compliance: kopen,k close .

Figure 3 .

 3 Figure 3.3 represents the general variation of the valve state in time. If one considers δP to be a constant in the set of equations (3.2), the trend and the solution of the ODE can be estimated.Then in order to obtain the valve orifice diameter (D v ), the following equation is used. Here, D [i-1] and D [i] refers to the diameters of the two lymphangions surrounding the valve.

Figure 3 . 4 :

 34 Figure 3.4: Sketch of a lymphangion, exposing the valve and the container part with the list of those variables related to the contraction model.

Figure 3 . 5 :

 35 Figure 3.5: Graph of the constant values for kc (i) and kr (ii) for different value of ecen and a. For these graphs a = emax m 0 and b = e min m 0 , with 1 > a > b > 0. Since b is not large and does vary a lot, it is not considered as a parameter.

Figure 3 . 6 :

 36 Figure 3.6: Evolution of the contraction intensity in function of time. The left graph is the contracting process. The right graph is the relaxing process. The timescale are not relevant as they are dependent on ω, and the contracting or relaxing compliance: kc,kr. The working conditions of the contraction strength are define between the bottom green line ( e min m 0 ), and the top green line ( emax m 0 )

Figure 3 . 7 :

 37 Figure 3.7: Passive pressure as a function of the diameter of the lymphangion. The range of diameter has been selected from minimum to maximum of the simulation used in this section.

Figure 3 . 8 :Figure 3 . 9 :

 3839 Figure3.8: Active pressure in function of the diameter of the lymphangion. The pressure here is in Barye and the diameter in centimeter. ξ is the contraction intensity. The range of diameter has been selected from minimum to maximum of the simulation completed here.

  Figure 3.12 Figure 3.13

Figure 3 . 12 :

 312 Figure3.12: Evolution of the valve of the second lymphangion in time. The graphs from top to bottom represent the pressure difference across the valve (i & ii), the valve state (iii) and the diameter of the valve orifice (iv). When the pressure difference becomes negative, the valve state switches from 0 (open) to 1 (close) which leads to the valve shutting closed.

F

  Figure 3.14: Evolution of the passive and active pressure of the middle lymphangion in time.Graphs from top to bottom represent pressure of the second lymphangion (i), diameter (ii). (iii) is a standard pressure-volume diagram, it allows to see the work and path taken by the contraction and relaxation cycle.

Figure 3 .

 3 Figure 3.15: Schematic of a single lymphatic vessel, made of two valves. This geometry represents the experimental set up made by Davis et al. [2012]. Here, PL represents the pressure inside the lymphangion. The rectangles define the position of the camera to record both the central lymphangion diameter as well as the n + 1.

  Maximum and minimum stiffness (e max & e min ) are mostly important in the contraction phase to define to which extend the contraction term of the constitutive equation can grow. Additionally, they define the overall period of a contraction-relaxation cycle since they define when to switch from one to another (cf. equation (3.4)). P d and D d control the stress-free condition for the constitutive equation. Both stress-free state pressure (P d ) and diameter (D d ) are used as parameters in equations (3.9) & (3.10). m 0 is the key constant in the active pressure relation as it sets the pressure scale. The couple of constants c 1 and c 2 are defined in the frequency equation and link the active pressure to diameter. They only matter when the pressure difference forces the fluid in the lymphangion and adapts the active response. f min is the minimum frequency, it sets the base level on the term ω of the active pressure equation (3.10).

  Figure 3.18: Graphs representing pressure (Ba) against time (s) for the simulated data (i)

  Figure 3.19: Schematic representation of the data set number 2 showing a lymphatic vessel in a chain of 5 lymphangions (i). Graph of the evolution of the boundary conditions (ii).

Figure 3 .

 3 Figure 3.20: Experimental data of a lymphatic vessel made out of five lymphangions. Plots are representing the diameter (cm) in function of time (s). From top to bottom, are the first to last lymphangions displayed.

Figure 3 .

 3 Figure 3.21: Simulation of a lymphatic vessel made out of five lymphangions. Plots are showing diameter(cm) in function of time(s). Top graphs represent the first lymphangion of the chain until the last one represented by the last line of graphs.

Figure 3 . 22 :

 322 Figure 3.22: Schematic of a simple network of three vessels, in a confluence configuration (left) and bifurcation configuration (right). Here the number of lymphangion per vessel (n) is three. However, this parameter varies for this study.

  Figure 3.23: Time average work per contraction for multiple n configurations, with n from three to ten. (i) graph represent the vc1 vessel while the (ii) graph the vc3 vessel.

Figure

  Figure3.23 (i) shows that the afferent vessels tend does not balance the work between each lymphangion properly. In fact, the lymphangions closer to the merging node produces more work than the one close to the input of the vessel. It is less the case in the (ii) graph, as one can observe that the work is approximately the same for all lymphangions from the confluence node to the output.

Figure 3 . 24 :

 324 Figure 3.24: Work per lymphangion against different boundary conditions for different length of vessels from 3 to 10 lymphangions. (i) graph represents the afferent top vessel (vc1), This vessel is efferent (vc3) on the (ii) graph.

Figure 3

 3 Figure3.24 (i) can be distinguished into two main groups: vessels with great number of lymphangions (n ≥ 7) and vessels with less than seven (n < 7). The first one grows more or less linearly with the pressure difference on the boundary conditions. However, the second group behaves differently. Indeed, for most cases when δP BC is set to 4 cm H 2 O, a maximum appears and configuration with three lymphangions displays a peak work. Also four lymphangions configuration peak work is shown for δP BC of 7 cm H 2 O.In figure3.24 (ii), the overall behaviour of this vessel is quite different from the previous one. It shows that whatever the number of lymphangion per vessel, the greater δP BC the larger work is provided in order to push lymph. It is possible to observe as well that for large δP BC (7 cm H 2 O and above) the behaviour is quite similar. A simple explanation could be that more lymphangions per vessel result in lower to individual pumping work. Besides, for large value of n and small values of δP BC , the average work appears to be constant, which is the base level the system remains when small difference of pressure are applied, such as observed in the previous section.Figure3.25, shows the average work per lymphangion multiplied per n versus the boundary condition pressure gradient. For the afferent vessel, a linear trend can be observed. With a pressure difference of 12 cm H 2 O, the average work per lymphangion multiplied per n grows by 0.01 erg approximately. The slope is just about 1.2 × 10 3 cm H 2 O erg -1 . But for figure 3.25 (ii), the curve resemble to a power law function or an exponential. This similarity is tested in for

Figure 3 .Figure 3 Figure 3 . 26 :

 33326 Figure3.24 (i) can be distinguished into two main groups: vessels with great number of lymphangions (n ≥ 7) and vessels with less than seven (n < 7). The first one grows more or less linearly with the pressure difference on the boundary conditions. However, the second group behaves differently. Indeed, for most cases when δP BC is set to 4 cm H 2 O, a maximum appears and configuration with three lymphangions displays a peak work. Also four lymphangions configuration peak work is shown for δP BC of 7 cm H 2 O.In figure3.24 (ii), the overall behaviour of this vessel is quite different from the previous one. It shows that whatever the number of lymphangion per vessel, the greater δP BC the larger work is provided in order to push lymph. It is possible to observe as well that for large δP BC (7 cm H 2 O and above) the behaviour is quite similar. A simple explanation could be that more lymphangions per vessel result in lower to individual pumping work. Besides, for large value of n and small values of δP BC , the average work appears to be constant, which is the base level the system remains when small difference of pressure are applied, such as observed in the previous section.Figure3.25, shows the average work per lymphangion multiplied per n versus the boundary condition pressure gradient. For the afferent vessel, a linear trend can be observed. With a pressure difference of 12 cm H 2 O, the average work per lymphangion multiplied per n grows by 0.01 erg approximately. The slope is just about 1.2 × 10 3 cm H 2 O erg -1 . But for figure 3.25 (ii), the curve resemble to a power law function or an exponential. This similarity is tested in for

Figure 3 . 27 :

 327 Figure 3.27: Time average work per contraction for multiple n configuration. (i) graph represent the v b1 vessel while the (ii) graph the v b3 vessel.

Figure 3 Figure 3

 33 Figure3.28: Work per lymphangion against different boundary conditions for different length of vessels from 3 to 10 lymphangions. On the left, vessel considered here is the afferent one (v b1 ), on the right vessel is the efferent one (v b2 ).

Figure 3

 3 Figure 3.30: Contraction per lymphangion per minute for different boundary conditions for different length of vessels from 3 to 11 lymphangions. Left graph represents the afferent vessel (v b1 ). The vessel represented on the right is the efferent one (v b2 ).

  Digitisation & general network structureFigures 3.31 & 3.32, show the original images and results of the digitisation. The process is as follows, the e cen e max e min D d

  Figure 3.31: (i) Example of lymphatic networks provided by figure 1.2. This drawing will be referenced as drawing one. Here blue dots represent bifurcation in the network and yellow represent confluences. The right red line represent the beginning of the network considered here, and left red line the end. The green lines can represent the inputs or outputs depending on their position. This is the first leg lymphatic network drawn in Sappey [1874].(ii) Disposition of the vessels replicating the above network with vessel numbers. The color green refers to the inlet, and red the outlets. Bifurcations and confluences are blue and yellow respectively.

  Figure 3.32: (i) Example of lymphatic networks provided by figure 1.2. This drawing will be referenced as drawing two. Here blue dots represent bifurcation in the network and yellow represent confluences. The right red line represent the beginning of the network considered here, and left red line the end. The green lines can represent the inputs or outputs depending on their position. This is the second leg lymphatic network drawn in Sappey [1874]. (ii) Disposition of the vessels replicating the above network with vessel numbers. The color green refers to the inlet, and red the outlets. Bifurcations and confluences are blue and yellow respectively.

Figure 3 . 33 :

 333 Figure 3.33: Pressure for all lymphangions against time, for the 61st vessel of the first drawn network.

Figure 3 . 34 :

 334 Figure 3.34: Pressure for all lymphangions against time, for the 17th vessel of the network from the first drawing.

Figure 3 . 35 :

 335 Figure 3.35: Pressure for all lymphangions against time, for the 103rd vessel of the second drawn network.

Figure 3 . 36 :

 336 Figure 3.36: Pressure for all lymphangions against time, for the 72nd vessel for the network of the second drawing.

Figure 3 . 37 :

 337 Figure 3.37: Pressure in the network in the initialisation phase, in the network of the drawing one presented in figure 3.31

Figure 3 . 38 :

 338 Figure 3.38: Pressure in the network in the initialisation phase, in the network of the drawing two presented in figure 3.32

Figure 4 . 1 :

 41 Figure 4.1: Flow chart of the main steps for the ISPM (Immersed Structural Potential Method) code. The top represents the initialisation and the bottom the end of the procedure. In green rhombi are the conditions occurring in the code logic. Here, t represents time, ∆t

Figure 4 . 2 :

 42 Figure 4.2: Staggered grid for fluid computation, representing a pressure centered computational cell. Here i & j are positive integers referring to the cell position. Pi,j represents the pressure and is estimated at the center of the cell i,j. The x-axis velocity, u is computed on the left and right cell's wall. Same for v the y-axis on the top and bottom.

Figure 4 . 3 :

 43 Figure 4.3: Representation of the Peskin [2002] kernel in function to its center. (i) Shows the kernel function φ in function of its radius. (ii) Displays the 2D representation of the kernel function.

Figure 4 . 4 :

 44 Figure 4.4: (i) Shear flow model used for numerical simulations of solid particle deforming in a Couette flow. With l the size of the simulation box, and r the radius of the particle simulated. (ii) Modification to the particle submitted to shear flow. Here a and b are the semi-major and semi-minor axis of the ellipse.

Figure 4 . 5 :

 45 Figure 4.5: Relative error of the velocity field and the pressure for a particle in a shear flow. The error for the velocity in the x-direction, y-direction as well as the pressure field are compared to the finest mesh case. This case is the run number 4 from table 4.1. Top and bottom slopes are here to represent respectively the error in O(h) and O(h 2 ).

Figure 4 . 6 :Figure 4 . 7 :

 4647 Figure 4.6: (i) -Representation of the particle solid kernels in a shear flow after a hundred timesteps, one can observe the beginning of the deformation. (ii) -Shape of the particle under shear flow after saturation of the deformation at time step t = 0.13 s only showing the solid kernels at the fluid-solid interface and not the whole set of solid kernels.

Figure 4 . 8 :

 48 Figure 4.8: Representation of the simulation setup in order to test the dispersion of particles in a Poiseuille flow.

Figure 4 .

 4 Figure 4.8 represents the geometry of the simulation. The goal of such a geometry is to evaluate the effect of the presence of the valve on a pulsed flow with particle. In order to achieve this, one needs to run multiple cases with different configurations. This batch of simulations presented here is made of seven specific setups. All these simulations use the dimensions defined in table4.2. When compared with the dimension of a real lymphangion only half a lymphangion length is simulated in order to save up computational time. It is a limit of this study to not have proper sized lymphangion, as it would contribute to dispersion of particle to a greater extent. Particles are chosen to fit the lymphocyte size reported in literature. As the number of particles varies with the immune response, it is interesting to test different particle concentrations. In this study, we focus on three cases with respectively 1, 36 or 72 particles, which corresponds to a density of approximately 2, 80 and 160 particles/mm 2 respectively. Furthermore, it is important to notice that the initial position of particles in the simulation box is random. For this geometry, a set of seven simulations is summarized in tables 4.3 and 4.2.In terms of dimensionless numbers the Reynolds, the Womersley and the Strouhal numbers are defined as:

Figure 4 .

 4 Figure 4.10: The error for the velocity in the x-direction, y-direction as well as the

Figure 4 .

 4 Figure 4.11 shows the entire valveless simulation domain at different capture times. The reader can appreciate the proportions as well as the size of particles and how they spread in the channel. One can observe that the particle tends to gather in the center and distribute themselves properly across all the entire channel. Also, for capture time 22.5 s, the Poiseuille flow is going on the opposite direction, since the input velocity frequency is 0.5 Hz as shown in figure4.9 as u in is a periodic function. This frequency is higher compare the contraction frequency showed in the lumped model chapter, because the simulation does

Figure 4 .

 4 Figure 4.11: Velocity fields of a pulsed flow with a particle representing a lymphocyte. Per column: (i) one particle simulation -case B, (ii) 36 particles -case C, (iii) 72 particlescase D. Per line: different capture times: 7.5, 15, 22.5 and 30 s.

Figure 4 .

 4 Figure 4.12: Velocity fields of a pulsed flow with a particle representing a lymphocyte and valve. Per column: (i) one particle simulation -case E, (ii) 36 particles -case F, (iii) 72 particles -case G. Per line: different capture times: 1.25, 2.5, 3.75 and 5 s.

Figure 4 . 13 :

 413 Figure 4.13: Distribution of particles in simulation without valve along the y-axis for different times for the simulation C, with 36 particles (i). And for the simulation D with 72 particles (ii), both without valve.

Figure 4 .Figure 4 . 14 :Figure 4 . 15 :

 4414415 Figure 4.14: Distribution of particles in simulation with valve along the y-axis for different times for the F simulation with 36 particles (i), and simulation G with 72 particles (ii) are present and are spread randomly at time 0 s.

Figure 4 . 16 :

 416 Figure 4.16: Distance of the closest particle from the top or bottom wall as a function of time for the simulation without (i) and with (ii) valves.

Figure

  Figure 4.16 shows the distance in the y-direction between the walls and its closest particle. This is a measure of the thickness of the particle-free layer near walls. (i) represents only valveless simulations, and (ii) simulations with valve. In the valveless cases, the observation formulated in the previous paragraph are

Figure 4 . 17 :

 417 Figure 4.17: Evolution of the position of y-coordinate particles in function of time. 0.015 cm is the middle of the channel and the Poiseuille flow is defined with y-orientation. In these figures, there are 36 particles (i) and 72 (ii) in the simulation so respectively C & D.

Figure 4 .

 4 Figure 4.18: Evolution of the position of a y-coordinate particles in function of time for valved simulation. In these figures, there are 36 particles (i) and 72 (ii) so respectively simulations F & G.

Figure 4 . 19 :

 419 Figure 4.19: Distribution of particle y-velocity (migration velocity) averaged in space and time for simulation without valve along the y-axis for different time ranges for the simulation C, with 36 particles (i). And for the simulation D with 72 particles (ii), both without valve.

Figure 4 .

 4 Figure 4.19, particles tends to get closer to the center (low shear rate zone) as

Figure 4 . 20 :

 420 Figure 4.20: Distribution of particle y-velocity (migration velocity) averaged in space and time for simulation with valve along the y-axis for different time ranges for the F simulation with 36 particles (i), and simulation G with 72 particles (ii). Particles in both cases are present and are spread randomly at initial time.

Figure 4 .Figure 4 .

 44 Figure 4.21: x-velocity profile as a function of y, for different times and different y-position. Each curve color represent a time in second given in the legend. (i) Figure uses results from the C simulation. (ii) Figure uses results from the D simulation.

Figure 4 .

 4 Figure 4.22: x-velocity profile as a function of y, for different times and different y-position. Each curve color represent a time in second given in the legend. (i) Figure uses results from the F simulation. (ii) Figure uses results from the G simulation.

Figure 4 . 23 :

 423 Figure 4.23: Evolution of the lower valve tip position. Here, the difference comes from the number of particles in the simulation. Two arrows represent the time evolution.

Figure 4 . 24 :

 424 Figure 4.24: Zoom on the solution of the lower valve tip position between 2 and 5 seconds of physical time. The difference comes from the number of particles in the simulation. The arrow shows the evolution of the curve as function of time.

Figure 4 . 25 :

 425 Figure 4.25: Representation of the simulation setup in order to understand the wall impact on lymph flow in a lymphangion.

Figure 4 .

 4 Figure 4.25 displays the geometry of this section simulations. The straight bars in the simulation box are walls, and are considered as solid, thus they interact with the fluid. The valves are identical to the previous section, and the position overlaps slightly with the top and bottom wall. So they are "attached" virtually, even if they are defined as different in the simulations.In table4.5, the values of the geometry used for this section 's simulations are summed up. The length is defined based on the half lymphangion geometry from the previous simulations. The height is defined to allow the movement of the wall on the y-direction. Leaflets and walls have been chosen based on the geometry presented by[START_REF] Daniel | Integrated geometric and mechanical analysis of an image-based lymphatic valve[END_REF] andZawieja [2009]. The geometry proposed here imitates the experimental setup presented in appendix. Its main objective is to move the walls as defined in figure4.25 to compress or expand the lymphatic volume. In order to achieve this goal, multiple strategies were used.

Figure 4 .

 4 Figure 4.26: Relative error of the velocity field and the pressure for the wall and valve study. The error for the velocity in the x-direction, y-direction as well as the pressure field are compared to the finest mesh case. This case is the run number 4 from table 4.6. Top and bottom slopes are here to represent respectively the error in O(h) and O(h 2 ).

Figure 4 .

 4 Figure 4.27: Initial geometry with wall pushing fluid outward. The solid structure is independent of the fluid domain, so it begins to be advected. Physical time is 1.07 s.

Figure 4 .

 4 Figure 4.28: Initial geometry with wall pushing fluid inward. The walls are submitted to forces pulling them apart from one another. Physical time is 0.55 s, since at 1 s this simulation diverge and the walls are submitted to instabilities.

Figure 4 .

 4 Figure 4.29: Second geometry with wall pushing fluid outward. In this case, the walls are attached to the simulation domain by constraining the x-velocity to zero where the wall and valve attaches.

Figure 4 .Figure 4 .

 44 Figure 4.30: (i) Second geometry flow at the left wall and the right wall. (ii) Sum of the previous flow in order to estimate the flow produced.

Figure A. 2 :

 2 Figure A.2: Dimension of a microfluidic collecting lymphangion like contractile body.

Figure

  Figure B.2 describes a basic lymphangion with relevant physical quantities. This figure sets up the physical problem to estimate the pressure profile of a

  integrating over the x-axis between r and -r, one should obtain: V a from the equation (B.22) into equation (B.24) one gets the following relation: r 3 η ∂ 2 p ∂y 2 + u = 0 (B.25)

  (i) Picture representing lymph and blood circulatory system and how they connect with each other. (ii) Inner leg lymphatic skin network from large drawing of Sappey[1874]. (iii) Microscopic image showing the profile section of a lymphangion with the two leaflets from[START_REF] Moore | Lymphatic system flows[END_REF]. . . . . . . . . . . . . . 1.2 Lymphatic vessels of the leg skin (Vaisseaux lymphatiques superficiels ou cutanés du membre inférieur in french) this illustration can be found in Sappey [1874]. . . . . . . . . . . . . . . . . . . 1.3 Schematic representation of one lymphangions with its neighboring ones. The Diamter here represent the lymphangion diameter, and the length shows the overall length of a lymphangion. The Valves are shown in two different state. They are oriented in a direction allowing it to close if the pressure difference moves the fluid backward, and close the valve otherwise if it moves forward. 1.4 Summary of elasticity values (Young's modulus) for different immune cells, from Bufi et al. [2015]. . . . . . . . . . . . . . . . . 2.1 Schematic of a network used in Jamalian et al. [2016] paper, here it is drawn with three lymphangions per vessels. In the article, the amount of lymphangions per vessel is varied. Furthermore, the overall structure remains the same, such that the number of connection between vessels remains at three with four inlets and one outlet. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Different types of valve studied from a numerical bio engineering point of view. (i) Aortic valve for numerical simulation from Spühler et al. [2018], (ii) idem from Loon, van [2005], (iii) Lymphatic valve from experimental result from Watson et al. [2017].3.1 Sketch of a lymphangion's side view, exposing the valve and the container part. (i) The position of different pressure points and diameters are provided. As well as the different flux between the pressure points (ii) The short arrows represent the wall motion of a contracting ith lymphangion. The short arrows in the i+1th lymphangion show the valve opening. Large arrows display the different pressures a lymphangion is submitted to. . . . . . . . . 3.2 Sketch of a lymphangion's valve, with the mesh pressure point. Table of the constant relevant for the valve model. Here, δP = P v[i] -P n[i-1] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.3 Evolution of the valve state in function of time. The left graph is the closing process. The right graph is the opening process. The time scales are not relevant as they are dependent on δP , and the opening or closing compliance: k open ,k close . . . . . . . 29 3.4 Sketch of a lymphangion, exposing the valve and the container part with the list of those variables related to the contraction model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.5 Graph of the constant values for k c (i) and k r (ii) for different value of e cen and a. For these graphs a = emax m0 and b = emin m0 , with 1 > a > b > 0. Since b is not large and does vary a lot, it is not considered as a parameter. . . . . . . . . . . . . . . . . . 31 3.6 Evolution of the contraction intensity in function of time. The left graph is the contracting process. The right graph is the relaxing process. The timescale are not relevant as they are dependent on ω, and the contracting or relaxing compliance: k c ,k r . The working conditions of the contraction strength are define between the bottom green line ( emin m0 ), and the top green line ( emax m0 ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.7 Passive pressure as a function of the diameter of the lymphangion. The range of diameter has been selected from minimum to maximum of the simulation used in this section. . . . . . . . . . 34 3.8 Active pressure in function of the diameter of the lymphangion. The pressure here is in Barye and the diameter in centimeter. ξ is the contraction intensity. The range of diameter has been selected from minimum to maximum of the simulation completed here. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.9 Sum of the active and the passive pressure in function of the diameter of the lymphangion. The pressure here is measured in Barye and the diameter in centimeter. The range of diameter has been selected from minimum to maximum of these simulations. 35 3.10 Zero-dimension approximation of a lymphangion using the electrical analogy. Dashed lines represent the limit of a lymphangion, diodes are the valve and capacitor represent the container/vessel. 36 3.11 Schematic of a single lymphatic vessel, made of three lymphangions. This set up is simple enough to test different parameters such as the valve or the contraction model. . . . . . . . . . . . . 42 3.12 Evolution of the valve of the second lymphangion in time. The graphs from top to bottom represent the pressure difference across the valve (i & ii), the valve state (iii) and the diameter of the valve orifice (iv). When the pressure difference becomes negative, the valve state switches from 0 (open) to 1 (close) which leads to the valve shutting closed. . . . . . . . . . . . . . . . . 43 3.13 Evolution of the contraction of the second lymphangion in time.

  The lymphangion pressure is shown to compare the contribution from the contraction pressure term. . . . . . . . . . . . . . . . 3.14 Evolution of the passive and active pressure of the middle lymphangion in time. Graphs from top to bottom represent pressure of the second lymphangion (i), diameter (ii). (iii) is a standard pressure-volume diagram, it allows to see the work and path taken by the contraction and relaxation cycle. . . . . . . . . . . 3.15 Schematic of a single lymphatic vessel, made of two valves. This geometry represents the experimental set up made by Davis et al. [2012]. Here, P L represents the pressure inside the lymphangion. The rectangles define the position of the camera to record both the central lymphangion diameter as well as the n + 1. . . . . 3.16 Schematic representation of the geometry for data set number 1 (i). Graph showing the evolution of pressure boundary conditions (ii). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.17 Diameter (cm) of second lymphangion vs time (s) for the simulation (i) and the experiment (ii). In terms of frequency, one can observe faster contractions when the pressure difference is larger and amplitude is reduced. Maximum of peak values are constant in the simulation while they decrease in the experimental data. 3.18 Graphs representing pressure (Ba) against time (s) for the simulated data (i) and the digitized data (ii). Here the envelop matches well between the two graphs. The frequency trend is properly captured as well. . . . . . . . . . . . . . . . . . . . . . 3.19 Schematic representation of the data set number 2 showing a lymphatic vessel in a chain of 5 lymphangions (i). Graph of the evolution of the boundary conditions (ii). . . . . . . . . . . . . . 3.20 Experimental data of a lymphatic vessel made out of five lymphangions. Plots are representing the diameter (cm) in function of time (s). From top to bottom, are the first to last lymphangions displayed. . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.21 Simulation of a lymphatic vessel made out of five lymphangions. Plots are showing diameter(cm) in function of time(s). Top graphs represent the first lymphangion of the chain until the last one represented by the last line of graphs. . . . . . . . . . . 3.22 Schematic of a simple network of three vessels, in a confluence configuration (left) and bifurcation configuration (right). Here the number of lymphangion per vessel (n) is three. However, this parameter varies for this study. . . . . . . . . . . . . . . . . 3.23 Time average work per contraction for multiple n configurations, with n from three to ten. (i) graph represent the v c1 vessel while the (ii) graph the v c3 vessel. . . . . . . . . . . . . . . . . . . . . 53 3.24 Work per lymphangion against different boundary conditions for different length of vessels from 3 to 10 lymphangions. (i) graph represents the afferent top vessel (v c1 ), This vessel is efferent (v c3 ) on the (ii) graph. . . . . . . . . . . . . . . . . . . . . . . . 54 3.25 Work per lymphangion multiply by the number of lymphangion in the vessel against different boundary conditions for different length of vessels from 3 to 10 lymphangions. (i) graph represents the afferent top vessel (v c1 ), This vessel is efferent (v c3 ) on the (ii) graph. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.26 Contraction per lymphangion per minute against different boundary conditions for different length of vessels from 3 to 11 lymphangions. (i) graph refers to vessel (v c1 ) and is the afferent bottom one, This vessel is the efferent one (v c3 ) on the (ii) graph. 55 3.27 Time average work per contraction for multiple n configuration. (i) graph represent the v b1 vessel while the (ii) graph the v b3 vessel. 56 3.28 Work per lymphangion against different boundary conditions for different length of vessels from 3 to 10 lymphangions. On the left, vessel considered here is the afferent one (v b1 ), on the right vessel is the efferent one (v b2 ). . . . . . . . . . . . . . . . . . . . 57 3.29 Work per lymphangion multiplied by the number of lymphangion in the vessel against different boundary conditions for different length of vessels from 3 to 10 lymphangions. (i) the vessel considered here is the afferent one (v b1 ), (ii) the vessel is the efferent one (v b2 ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.30 Contraction per lymphangion per minute for different boundary conditions for different length of vessels from 3 to 11 lymphangions. Left graph represents the afferent vessel (v b1 ). The vessel represented on the right is the efferent one (v b2 ). . . . . . . . . 58 3.31 (i) Example of lymphatic networks provided by figure 1.2. This drawing will be referenced as drawing one. Here blue dots represent bifurcation in the network and yellow represent confluences. The right red line represent the beginning of the network considered here, and left red line the end. The green lines can represent the inputs or outputs depending on their position. This is the first leg lymphatic network drawn in Sappey [1874]. (ii) Disposition of the vessels replicating the above network with vessel numbers. The color green refers to the inlet, and red the outlets. Bifurcations and confluences are blue and yellow respectively. . 60 3.32 (i) Example of lymphatic networks provided by figure 1.2. This drawing will be referenced as drawing two. Here blue dots represent bifurcation in the network and yellow represent confluences. The right red line represent the beginning of the network considered here, and left red line the end. The green lines can represent the inputs or outputs depending on their position. This is the second leg lymphatic network drawn in Sappey [1874]. (ii) Disposition of the vessels replicating the above network with vessel numbers. The color green refers to the inlet, and red the outlets. Bifurcations and confluences are blue and yellow respectively. . 3.33 Pressure for all lymphangions against time, for the 61st vessel of the first drawn network. . . . . . . . . . . . . . . . . . . . . . . 3.34 Pressure for all lymphangions against time, for the 17th vessel of the network from the first drawing. . . . . . . . . . . . . . . 3.35 Pressure for all lymphangions against time, for the 103rd vessel of the second drawn network. . . . . . . . . . . . . . . . . . . . 3.36 Pressure for all lymphangions against time, for the 72nd vessel for the network of the second drawing. . . . . . . . . . . . . . . 3.37 Pressure in the network in the initialisation phase, in the network of the drawing one presented in figure 3.31 . . . . . . . . . . . . 3.38 Pressure in the network in the initialisation phase, in the network of the drawing two presented in figure 3.32 . . . . . . . . . . . . 4.1 Flow chart of the main steps for the ISPM (Immersed Structural Potential Method) code. The top represents the initialisation and the bottom the end of the procedure. In green rhombi are the conditions occurring in the code logic. Here, t represents time, ∆t is the time-step, Err(f FSI) = |f n FSI -f n-1 FSI | |f n FSI |the error function and k tolerance is the tolerance constant. . . . . . . . . . . . . . 4.2 Staggered grid for fluid computation, representing a pressure centered computational cell. Here i & j are positive integers referring to the cell position. P i,j represents the pressure and is estimated at the center of the cell i,j. The x-axis velocity, u is computed on the left and right cell's wall. Same for v the y-axis on the top and bottom. . . . . . . . . . . . . . . . . . . . . . . 4.3 Representation of the Peskin [2002] kernel in function to its center. (i) Shows the kernel function φ in function of its radius. (ii) Displays the 2D representation of the kernel function. . . . . . 4.4 (i) Shear flow model used for numerical simulations of solid particle deforming in a Couette flow. With l the size of the simulation box, and r the radius of the particle simulated. (ii) Modification to the particle submitted to shear flow. Here a and b are the semi-major and semi-minor axis of the ellipse. . . . . . . . . .4.26 Relative error of the velocity field and the pressure for the wall and valve study. The error for the velocity in the x-direction, y-direction as well as the pressure field are compared to the finest mesh case. This case is the run number 4 from table 4.6. Top and bottom slopes are here to represent respectively the error in O(h) and O(h 2 ). . . . . . . . . . . . . . . . . . . . . . . . . . . 101 4.27 Initial geometry with wall pushing fluid outward. The solid structure is independent of the fluid domain, so it begins to be advected. Physical time is 1.07 s. . . . . . . . . . . . . . . . . . 102 4.28 Initial geometry with wall pushing fluid inward. The walls are submitted to forces pulling them apart from one another. Physical time is 0.55 s, since at 1 s this simulation diverge and the walls are submitted to instabilities. . . . . . . . . . . . . . . . . 102 4.29 Second geometry with wall pushing fluid outward. In this case, the walls are attached to the simulation domain by constraining the x-velocity to zero where the wall and valve attaches. . . . . 103 4.30 (i) Second geometry flow at the left wall and the right wall. (ii) Sum of the previous flow in order to estimate the flow produced.104 4.31 (ii) Diameter here represent the distance between the top and bottom wall, evaluated close to the wall. . . . . . . . . . . . . . 104 A.1 Dimension of a microfluidic collecting lymphangion like valve. . 111 A.2 Dimension of a microfluidic collecting lymphangion like contractile body. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 B.1 Schematic of a lymphangion's without taking the valve into account. Here, L represents the length of a lympahngion, u the velocity field of the lymph in a Poiseuille flow, h the diameter of the lymphangion. . . . . . . . . . . . . . . . . . . . . . . . . . . 115 B.2 Schematic of a lymphangion's chain with a symmetry showing the upper half of the chain. Here, L represents the length of a lympahangion, u the velocity of contraction of the lymphangion upper wall, r the radius of the lymphangion and V a the mean velocity of the lymph in the lymphangion. . . . . . . . . . . . . 118

  List of different sources for dimensions of typical human lymphangion of the collecting network. (x) denotes that the value could not be retrieved or it is a variable.

	Margaris	Rahbar and	Jamalian	Bertram
	and	Black	Moore [2011]	et al. [2016]	[2020]
	[2012] Pan		Contarino
	et al. [2010]		and	Toro
				[2018]
	Diameter 0.2 mm		0.2 mm	x		x
	Length 2 mm		2 mm	3 mm		x
	Wall thickness x		10 µm	x		x
	Valve thickness x		x	x		5 µm
	Table 1.1:

Table 1 .

 1 2: Characteristic values of a lymphocyte that affects lymph flow.

	Kristenson	Petchakup	Bufi et al.
	[2009]	et al. [2021]	[2015]
	Diameter 5 ∼ 7 µm	8 ∼ 10 µm	x
	Rigidity x	x	90 Pa

Table 2

 2

.1:

  Table of main different constitutive equations found in literature.

		Method	Geometry	Metrics
	Bertram et al. [2011]	Forced pulsation to con-trol contraction cycle Poiseuille flow assumption Working valve with differ-ent opening and closing Unknown variables are pressure difference	ODE chain of two valves and one lymphangion Forced contraction cycles.	• Pressure • Diameter • Valve resistance • Flow
		pressure and diameter.
	Contarino and Toro [2018]	Poiseuille flow assumption Valve opens and closes with respect to the pres-sure difference Unknown variables are area and flow.	PDE series of three lymphan-gions mechanico-chemical bal-ance for contraction cy-cles.	• Pressure • Diameter • Valve state • Flow • Stimulus, contrac-
				tion state
	Jamalian et al. [2013]	Poiseuille flow assumption for flow resistance Forced pulsation to con-trol contraction cycle Sensitivity analysis for rameters valve and contraction pa-	ODE Variable chain from 3 to 21 lymphangions Forced contraction cycles.	• Pressure • Diameter • Frequency • Flow
		Unknown variables are
		pressure and diameter.
	Tretyakova et al. [2018]	Poiseuille flow assumption for flow resistance Forced pulsation to con-trol contraction cycle Unknown variables are section and velocity Constant Poiseuille resis-	ODE and PDE chain of 2 to 10 lymphan-gions Lymphatic channel valves are treated differently from 0D to 1D.	• Pressure • Diameter • Flow • Active contraction term
		tance for valves in 0D
		Larger friction force to
		prevent back flow in 1D.

Table 2 .

 2 2: Recapitulating table of method and geometry from relevant literature papers on the secondary valve and lymphangion model.

Table 2

 2

	.3: Recapitulating table of method and geometry from relevant literature papers on
	the secondary valve and lymphangion 3D/2D model.

Chapter 3 Discrete modeling of lymphangion chains and networks

	Résumé -On présente un modèle d'écoulement dans le réseau lymphatique
	basé sur une modélisation 0D de la dynamique des lymphangion avec une nou-
	velle formulation numérique, pour le calcul d'écoulement de la lymphe dans le
	réseau collecteur. Pour cette étude, on se base sur une version simplifiée du
	modèle qui réduit le nombre de paramètres dans l'équation constitutive sans
	modifier sa physique. Cette première se portera sur différents cas spécifiques
	tel que la bifurcation, convergence, nécessaire pour la création d'un réseau. Les
	résultats seront comparés avec des données expérimentales effectués sur des
	rats de laboratoire. Enfin, on présente les résultats d'une géométrie de réseau
	spécifique extraite des planches d'anatomie d'une jambe humaine.

  Table of the constant relevant for the valve model. Here, δP = P v[i] -P n[i-1] .

  Table listing all variables and parameters for the constitutive equation.

	Passive Pressure
	P d pressure constant
	D d diameter constant
	Active Pressure
	τ thickness of the lymphangion wall (cm)
	D diameter (cm)
	γ contraction memory
	Table 3.2:

Table 3 .

 3 3: List of all parameters used for the 0d model, and their default or initial values.

	Parameter Description	Value	Units
		Boundary & Initial conditions
	P in	input pressure	3920	Ba
	P out	output pressure	4900	Ba
	P ext	external pressure	0	Ba
	P n[i,0]	initial pressure node	1960	Ba
	P v[i,0]	initial pressure valve	1960	Ba
		Contraction & Valve Mobilisation
	e cen	eccentricity of contraction	0.2	erg cm -1
	e min	min stiffness for the contraction function	10	erg cm -1
	e max	max stiffness for the contraction function	6e5	erg cm -1
	ω min	minimum of frequency	0.01	Hz
	m 0	stiffness factor	1e6	erg cm -1
	D d	contraction lymphangions diameter	2e-3	cm
	P d	contraction pressure for the constitutive law	35	Ba
	c 1	threshold for D/Dd to increase frequency	2.1	-
	c 2	rate of increase of frequency	0.2	-
	c	constant	2.15	-
	k open	valve opening compliance	100	Ba -1 s -1
	k close	valve closing compliance	100	Ba -1 s -1
		Material Mobilisation
	l	numbers of lymphangions -1	9	-
	L valve	valve length	0.01	cm
	ρ	lymph density	1	g cm -3
	µ	lymph dynamic viscosity	0.01	P
	D	lymphangions diameter	0.01	cm
	L	lymphangions length	0.3	cm
		Numeric Mobilisation
	δt	time step	5e-3	s
	T	simulation total time	10	s
	ε	epsilon	1e-10	-

  .4 defines bounding values used in this study. Eccentricity (e cen ) refers to the ratio of time spent between the contraction andTable 3.4: List of minimum and maximum values used to tune the contraction model against Davis et al. [2012] experimental data.

		Units	Minimum	Maximum
	eccentricity (e cen )	-	0.090	0.11
	max stiffness (e max )	-	2.7 × 10 5	3.3 × 10 5
	min stiffness (e min )	-	9.0 × 10 -5	1.1 × 10 -4
	stress-free contraction diameter (D d )	cm	2.5 × 10 -3	3.1 × 10 -3
	stress-free contraction pressure (P d )	Ba	26	32
	dimension active pressure (m 0 )	-	9.0 × 10 7	1.1 × 10 8
	zero stress active contraction (c)	-	1.71	2.09
	threshold to increase frequency (c 1 )	-	1.89	2.31
	rate of for frequency increase (c 2 )	-	4.5 × 10 -2	5.5 × 10 -2
	minimum frequency (f min )	Hz	9.0 × 10 -2	1.1 × 10 -1

Table 3

 3

.5: Values used to tune the model to fit the data set 1. Data set 2 values are shown in table 3.6.

Table 3 .

 3 6: Values used to tune the model to fit the data set 2. Data set 1 values are in table 3.5.

Table 3 . 7 :

 37 Values of the contraction constants used for the confluence and bifurcation network study.

Table 3 .

 3 8: Values used to model lymphangions in the full network simulations. Table of number of inputs, outputs, confluences and bifurcations for the two networks presented in both figures 3.31 & 3.32

		input output confluence bifurcation
	fig I.	72	25	62	15
	fig II.	46	26	39	16
	Table 3.9:

Table 3 .

 3 10: Details of the vessel showing chaotic behaviour studied in this section.

	Vessel	lymphangion	input	input type	output	output type
	number	number	position		position
	drawing 1	15	(32.42; 18.73) simulation	(34.38; 16.69)	confluence
	vessel 61			inlet

Table 3 .

 3 11: Details of vessels displaying in phase behaviour studied in this section.

	Vessel	lymphangion	input	input type	output	output type
	number	number	position		position
	drawing 1	16	(7.02; 6.27)	simulation	(10.08; 5.34)	confluence
	vessel 17			inlet
	drawing 2	19	(78.73; 19.56) simulation	(81.74; 17.49)	confluence
	vessel 103			inlet

Table 4 . 2 :

 42 List of the main physical values for the test of dispersion of particles in a Poiseuille flow.

	Simulation Domain input & output (l) 0.3 mm
					wall (L) 1.5 mm
			Leaflet		width 0.03 mm
					length ≈ 0.33 mm
	Boundary conditions			u in 1 mm s -1
			Particles		amount 1, 36 or 72
					diameter 0.01 mm
	Simulation	A	B	C	D	E	F	G
	Valve	with without without without with with with
	Number of	0	1	36	72	1	36	72
	particles

Table 4 .

 4 3: List of simulations used for this section geometry with their characteristics.

  It is described in table 4.4, run 4. Here the

	Run Timesteps/s	n x ,n y	Solid kernels	δt	h
	1	4.28 × 10 -2 69,21	48600	1.041 7×10 -6 2.173×10 -3
	2	4.06 × 10 -2 141,45	48600	5.208 3×10 -7 1.063×10 -3
	3	4.11 × 10 -2 285,93	48600	2.604 2×10 -7 5.263×10 -4
	4	4.03 × 10 -2 573,189	48600	1.302 1×10 -7 2.617×10 -4

  direction, y-direction as well as the pressure field are compared to the finest mesh case. This case is the run number 4 from table 4.6. Top and bottom slopes are here to represent respectively the error in O(h) and O(h 2 ).

	Run Timesteps/s	n x ,n y	Solid kernels	δt	h
	1	9.21 × 10 -3 31,21	135750	2.082 3×10 -5 4.838×10 -3
	2	9.37 × 10 -3 73,41	135750	1.044 7×10 -5 2.054×10 -3
	3	1.10 × 10 -2 145,81	135750	5.219 0×10 -6 1.034×10 -3
	4	1.13 × 10 -2 289,161	135750	2.603 6×10 -6 5.190×10 -4

Table 4 . 6 :

 46 List of the discretisation and performance parameters for the wall contraction and valve study. Here time-step per second refers to the amount of time-steps achieved in one seconds. δt correspond to the time discretisation of the simulation. nx and ny represent the number of fluid computational cells in each direction x and y. h is the length of a fluid cell, it is obtained by l nx and is mesured in centimeter.

  .1 shows the data set for the first leg lymphatic network represented in Sappey[1874]. And, tableC.2 shows the data set for the third leg lymphatic network represented in Sappey[1874].Table C.2:List of all nodes used to build the network from the third drawing representing the secondary leg lymphatic network.

	Tables
			13.6698	110 111
	3	88.5201	19.4744	112 114
	Outlet
	5	77.2337	1.27907
	5	77.1784	2.11628
	5	77.1231	2.84186
	5	77.2061	3.59535
	5	77.3444	4.48837
	5	77.7593	5.66047
	5	78.4786	7.11163
	5	80.3596	8.95349
	5	83.1535	10.0698
	5	85.2282	10.9349
	5	87.9391	11.5488
	5	89.8202	12.107
	5	91.2586	12.4977
	5	94.6058	13.4186
	5	24.7303	5.10233
	5	27.9115	4.51628
	5	32.3375	3.56744
	5	34.6611	3.0093
	5	42.296	0.69302
	5	44.9793	0.58139

  1.1 List of different sources for dimensions of typical human lymphangion of the collecting network. (x) denotes that the value could not be retrieved or it is a variable. . . . . . . . . . . . . . 1.2 Characteristic values of a lymphocyte that affects lymph flow. .2.1Table of main different constitutive equations found in literature. 2.2 Recapitulating table of method and geometry from relevant literature papers on the secondary valve and lymphangion model. 2.3 Recapitulating table of method and geometry from relevant literature papers on the secondary valve and lymphangion 3D/2D model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Table listing the units from the centimeter, gram and second system used in this document. . . . . . . . . . . . . . . . . . . . 3.2 Table listing all variables and parameters for the constitutive equation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 List of all parameters used for the 0d model, and their default or initial values. . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 List of minimum and maximum values used to tune the contraction model against Davis et al. [2012] experimental data. . . . . 3.5 Values used to tune the model to fit the data set 1. Data set 2 values are shown in table 3.6. . . . . . . . . . . . . . . . . . . . 3.6 Values used to tune the model to fit the data set 2. Data set 1 values are in table 3.5. . . . . . . . . . . . . . . . . . . . . . . 3.7 Values of the contraction constants used for the confluence and bifurcation network study. . . . . . . . . . . . . . . . . . . . . . 3.8 Values used to model lymphangions in the full network simulations. 3.9 Table of number of inputs, outputs, confluences and bifurcations for the two networks presented in both figures 3.31 & 3.32 . . . 3.10 Details of the vessel showing chaotic behaviour studied in this section. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.11 Details of vessels displaying in phase behaviour studied in this section. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.12 Details of the vessel showing appearance of synchronisation and studied in this section. . . . . . . . . . . . . . . . . . . . . . . .

NO is a relatively labile free radical with an in vivo half-life of less than 5 s, NO cannot be stored in free form and must generally be synthesized on demand or from more stable adducts (such as Glutathione) that have specific biological effects[START_REF] Dudzinski | The regulation and pharmacology of endothelial nitric oxide synthase[END_REF].

Using Dell ® Inspiron 3847 -Core i3 4150 3.5 GHz -8 GB -4 TB

Using Intel ® Xeon ® Gold 6148 CPU 2.4 GHz -80 CPU -100GB

This software is published by[START_REF] Mitchell | markummitchell/engauge-digitizer: Nonrelease[END_REF]

also known as kernel from[START_REF] Peskin | The immersed boundary method[END_REF]

Also known as loop unwinding, it is a computer science technique that a compiler attempts to optimize a program execution speed but enlarging the size of the binary file.

Their article use below 0.3 Ca as a limit.

Acknowledgements

This thesis has been produced from my work during three years, half spend in Swansea University and the other half in Université Grenoble Alpes. I was part of the team Biomedical Engineering research group in the ZCCE (Zienkiewicz Centre for Computational Engineering) and part of the MC2 (

Appendix A

Design of microfluidic experiments

Introduction This document aims at presenting microfluidic experiences. They can be used to understand flow at the scale of a collecting lymphangion. There are two different dynamic behaviours which are interesting from a modeling point of view. The first one is the valve, the second one is the vascular wall hereafter called contractile body. They will be both presented here, with dimensions and a rough estimation of the movement. Here, the microfluidic channel are made of PDMS gel on a microscope glass slide. The study of lymphangion using microfluidic channel is recent, [START_REF] Selahi | Abstract 274: 3d cylindrical lymphangion-on-a-chip: a new method to model lymphatic inflammatory & therapeutic responses[END_REF]. The approch taken in this cited abstract is more biological than the work presented here. However, the use of PDMS 1 -Glass microfluidic devices for micro-fluidic experiments is not recent. It is convinient to model and form into a micro channel structure. At the same time, it is biologically compatible and widely distributed. Using the static beam theory, one obtains the maximum deformation of a beam fixed at both ends with a uniform distribution using the following relation :

Where constant values are given in the table A.1.

uniform load ω 10 Pa to 50 Pa length L 0.3 mm modulus of elasticity E 1.5 MPa area moment of inertia I 9.5 × 10 -8 mm 3 The pressure difference between lymphangions, across a valve is in the range of 1 mm of water up to 5 mm, which is approximately 10 Pa to 50 Pa.

The formula to get the area moment of inertia is:

Using the geometry we have b = 30 µm and h ≈ 0.335 mm, thus I ≈ 9.5 × 10 -8 mm 3 . So we get the maximum deformation on a leaflet of the valve depending on the load of:

These results look quite small in regard to the size of the channel. For more interesting experience, it would be better to have a range ten times larger. This can be obtained by using a smaller modulus of 1.5 × 10 2 kPa or by thinning the width of the valve.

Contractile body

Here, the same approach for the valve and parameters is used but with a different I and ω. In fact, the pressure difference for the contractile walls should be about 1 bar, thus one obtains:

uniform load ω 1 × 10 5 Pa area moment of inertia I 8.2 × 10 -5 mm 3 According to figure A.2, we set the following parameters, b = 2.7 mm and h = 50 µm. So, the maximum displacement of the wall is:

This result gives a good approximation, but for experimental purposes, it would be better to have the estimation with an order of magnitude larger. A value closer to the half width of the channel would be better, which is around 0.15 mm.

Appendix B

Analytical solutions for a simplified lymphangion in 2D

Analytical solutions used to fit the solution in the mesh convergence study are defined here. A newtonian fluid is considered for all the following calculations.

Here, u = u x e x and defines the velocity vector field. p represents the pressure. ρ and µ are respectively the fluid density and the fluid viscosity. The bold font is used to represent any vector quantity.

The first equation represents the mass conservation of the fluid, and the second one the conservation of momentum. Often, a third equation is added for tracking a scalar field such as a concentration marker, or energy if temperature needs to be evaluated, which is the case of compressible fluid.

Velocity -Poiseuille flow

First, the Navier-Stokes equations (B.1) are written as a start ground for the derivation of the Poiseuille profile. The following hypothesis are considered:

• incompressible,

• newtonien fluid,

115

• stationnary flow,

• unidirectional flow.

The first one allows to consider ρ as a constant thus takes it out of the derivative. The second hypothesis simplifies the diffusion term of the Navierstokes momentum equation into

which can later be simplified into a laplacian operator. The third hypothesis removes any time-dependent derivative as the flow is considered established and steady thus time derivative is zero. The last hypothesis allows us to write the velocity only in one direction. Considering the problem geometry we will use it such as:

Also since the model for the lymphangion is in two-dimension, equations will also be developped in a two-dimension cartesian space. So by applying all the hypothesis, the mass conservation can be written:

Thus, one can deduce that the velocity only depends on the y-variable:

By using the conservation of momentum equation (B.1) in the x-direction, and by introducing the result of equation (B.5), it is possible to write the following equation:

By doing the same for the y-direction:

So the pressure only depends on x and not on y. One can then consider the hypothesis on the pressure. If the pressure is constant or does not depend on x, then it is possible to integrate equation (B.6) twice on the y-direction:

Here C 1 and C 2 are integration constants. They will be set using boundary conditions of the problem, such as detailled on figure B.1. Since the walls are considered to be no slip, one can write:

Then, one can set the values for the integration constant as:

So, the analytical solution for the velocity can then be written as:

Table of vessels mapping for full network zero dimension simulations

Introduction There are two drawings for secondary leg network provided by Sappey [1874]. Here, is the detail of all the vessels connection and number used in this thesis simulations. For simplicity purposes, foot is not taken into account in both drawing, so the network in this thesis start at the ankle and stops before the groin lymph nodes. Lymph vessels coming from the foot are referred to as lower inlets. Also conglomeration of capillaries lead to secondary vessels in the network, these vessels are called capillaries inlets Then, as the drawing only show the inner or outer face of different legs, some vessels can cross these two faces and they will be referred to as lateral inlets. With regard to outlets, they are only located before the groin area. In order to limit the amount of information process by the input file, only the inlet, outlet, bifurcation and confluence are defined, they will be referred to as connections. To sum up there are five different types of connections :

• Inlets are the entry points to the secondary lymphatic network and 1 is the connection number.

-Lower inlets from the secondary lymph vessels of the foot, -Capillary inlets from accumulation of capillaries, -Lateral inlets from the vessels coming from the other side of the leg.

• Bifurcations correspond a state where a vessel split into two vessels, its connection number is 2.

• Confluences represent a merging of two vessels into one, its connection number is 3.

• Links refer to a points linking a vessel to another, its connection number is 4, but not used here.

• Outputs correspond to the exit points of the network, 5 is the connection number.