Gregory Arbia

Chiara Corsini

Catriona Baker

Giancarlo Pennati

Tain-Yen Hsia

Irene E Vignon-Clementel
email: irene.vignon-clementel@inria.fr

Multiscale modeling of pulmonary hemodynamics in first stage single ventricle patients: integrating clinical data and simulations

Keywords: Parameter estimation, Clinical hemodynamics data, Patient-specific CFD, Multiscale simulations, First stage single ventricle patients, Pulmonary circulation modelling

Purpose: Congenital heart defects involve pathologies in which the heart has only one functional pumping chamber. In these conditions, treatments consist of three procedures where stage 1 occurs just after birth and links the systemic and pulmonary circulations through an artificial shunt. Computational fluid dynamics is a tool to explore hemodynamics in patientspecific simulations, and understand specific physiopathologies. Representative hemodynamics models are built based on imaging data for example to explore several virtual surgical options. In this context, the specification of boundary conditions is necessary to compute pressure and flow in the entire domain. However, these boundary conditions are rarely the measured variables. Besides, to take into account the rest of the circulation outside of the three-dimensional modeled domain, a number of reduced order models exist. Methods: In this work a method is presented to iteratively but automatically tune from clinically measured hemodynamics data reduced model parameters. Target pressures such as minimum, maximum or average over one or both lungs are considered, depending on the clinical measurement. A 5-parameter model representing more specifically the arterial and venous pulmonary circulation is proposed.

Results: The methodology is demonstrated with six stage 1 patient specific models: the clinical target data are well matched, except for one case where clinical data were incoherent. Inhomogeneous pressure, swirling blood flow patterns and very high wall shear stress 3D maps highlight similarities and differences among patients. Steady and pulsatile tuning results are similar. Conclusions: The highly complex fluid dynamics of this stage is for the first time analyzed around the distal systemic-to-pulmonary shunt anastomosis and the connected pulmonary ar-

Introduction

Single ventricle defects are a spectrum of congenital heart disorders resulting in only one functional ventricle. Several stages of palliative surgery are required to connect the systemic and pulmonary circulations to the single ventricular power source. The first stage commonly utilises a shunt between the aorta, or its branches, and the pulmonary arteries. Accurate modeling of patient-specific physiology for clinical decision-making requires the integration of the patient's clinical data into numerical simulations [START_REF] Vignon-Clementel | A primer on computational simulation in congenital heart disease for the clinician[END_REF]. Computational fluid dynamics studies of single-ventricle conditions [START_REF] Bove | Computational fluid dynamics in the evaluation of hemodynamic performance of cavopulmonary connections after the Norwood procedure for hypoplastic left heart syndrome[END_REF][START_REF] Vignon-Clementel | Outflow boundary conditions for 3D simulations of non-periodic blood flow and pressure fields in deformable arteries[END_REF][START_REF] De Zélicourt | Imaging and patient-specific simulations for the Fontan surgery: Current methodologies and clinical applications[END_REF] have presented the crucial interest of specifying the relevant boundary conditions. More specifically in [START_REF] Migliavacca | Modeling of the Norwood circulation: effects of shunt size, vascular resistances, and heart rate[END_REF] authors investigated the effect of clinical parameters such as shunt, vasculature and heart rate on Norwood procedure. A common way to take into account the effect of the downstream domain is to couple 3D Navier-Stokes equations to electric-analog (also called zero-dimensional or reduced) models or one-dimensional models of the peripheral circulation [START_REF] Pennati | Boundary conditions of patient-specific fluid dynamics modelling of cavopulmonary connections: possible adaptation of pulmonary resistances results in a critical issue for a virtual surgical planning[END_REF][START_REF] Vignon-Clementel | Outflow boundary conditions for three-dimensional finite element modeling of blood flow and pressure in arteries[END_REF][START_REF] Vignon-Clementel | Outflow boundary conditions for 3D simulations of non-periodic blood flow and pressure fields in deformable arteries[END_REF][START_REF] Vignon-Clementel | A primer on computational simulation in congenital heart disease for the clinician[END_REF]. This gives information about local fluid dynamics due to changes in anatomical features [START_REF] Corsini | An integrated approach to patient-specific predictive modeling for single ventricle heart palliation. Comput[END_REF][START_REF] Kung | Others: Predictive modeling of the virtual Hemi-Fontan operation for second stage single ventricle palliation: two patient-specific cases[END_REF]. In [START_REF] Qian | Computational Hemodynamic Analysis in Congenital Heart Disease: Simulation of the Norwood Procedure[END_REF] the authors presented computational hemodynamics simulations but in the systemic circulation of one specific Norwood case with aortic arch repair and Damus-Kay-Stansel anastomosis. The distribution of flows among the outlets was not directly controlled since a zero pressure gradient method was applied at the end of very long extensions to recover a physiological pressure level. The pressure results were compared with pressure measurements at three locations in the aorta.

A first idea consists of manually tuning parameters, which requires hemodynamics intuition [START_REF] Pennati | Boundary conditions of patient-specific fluid dynamics modelling of cavopulmonary connections: possible adaptation of pulmonary resistances results in a critical issue for a virtual surgical planning[END_REF] and has been experienced to be out of reach if there are many branches. Automatic parameter estimation methods have thus been developed. In [START_REF] Spilker | Tuning multidomain hemodynamic simulations to match physiological measurements[END_REF], six parameters of an aortic model were estimated with a quasi-Newton method to achieve some pressure and flow waveform features inside the 3D model. In [START_REF] Ismail | Adjoint-based inverse analysis of windkessel parameters for patient-specific vascular models[END_REF], the Windkessel parameters of 3D bifurcating aneurysm were estimated with an adjoint-based method to match systolic, diastolic and average pressure differences. In [START_REF] Bertoglio | Sequential parameter estimation for fluid-structure problems: Application to hemodynamics[END_REF], wall displacement values were used to estimate the stiffness of the 3D fluid-solid interaction idealized aneurysm model and its outlet proximal resistance, based on a sequential estimation approach. In [START_REF] Pant | A multiscale filtering-based parameter estimation method for patient-specific coarctation simulations in rest and exercise[END_REF], a dozen of Windkessel parameters in a 3D coarctation model were estimated with a Kalman filter approach on a corresponding 0D surrogate model of the 3D-0D model, to match given flow and pressure waveforms inside the 3D domain.

In this paper we present a CFD analysis of the blood flow behavior in the distal anastomosis of the systemic-to-pulmonary shunt and the connected pulmonary arteries over several bifurcations, in six different patients. At the 3D outlets, it would not be relevant to apply pressure boundary conditions for the following reasons: 1) the measured flow distribution between the two different lungs cannot be easily matched, and 2) the pressure measurements were not performed at these locations in the pulmonary arteries. Sometimes, applying the same pressure at all outlets can even lead to unphysiological reverse average flow in some of the branches. Applying timevarying flow boundary conditions would also be difficult because the flow distribution between the different outlets is not necessarily constant, and would not guaranty to be coherent with the pressure measurements. Moreover, prescribing pressure or flow reduces the predictive potential such as in virtual surgery planning. We thus aim at coupling 3D Navier-Stokes equations to reduced models including resistances and capacitances. A previous work [START_REF] Troianowski | Others: Three-dimensional simulations in Glenn patients: clinically based boundary conditions, hemodynamic results and sensitivity to input data[END_REF] introduced a method to iteratively tune reduced order parameters of the 3D model to match clinical inlet average pressure, inlet average flow and outlets' flow repartition. The method was illustrated on coupling of 3D steady Navier-Stokes equations to total outlet resistances of dozens of branches in Stage 2 (Glenn) patient-specific simulations. From these resistances, 3-parameter Windkessel models were generated to represent the downstream arterial tree, neglecting the venous circulation.

In this work we extend the method to iteratively tune reduced parameters for more general pressure clinical data localization, as clinical data are rarely given at the specific 3D boundaries. The methodology is demonstrated for steady and pulsatile cases, with clinical data of stage 1 patients. The highly complex fluid dynamics of this stage is analyzed. Moreover, the downstream reduced models of the pulmonary circulation is enriched to a 5-parameter model representing more specifically the arterial and venous circulations.

The article continues with the Methods section to present the different steps of the algorithm, from clinical data, automatic parameters tuning to the final reduced pulmonary model generation for each outlet. The Results section presents then the demonstration of this framework on six stage 1 single ventricle patients, reporting their specific 3D hemodynamics. A discussion follows on methodological, biomechanical and clinical aspects.

Methods

The Methods section presents the type of clinical data that define the patient-specific geometry and hemodynamics target values. The automatic reduced model parameters tuning to reach these target values is explained, followed by the 5-parameter reduced model description. In this work six patient-specific computational models were built to describe hemodynamics into the corresponding three-dimensional regions of interest, reflecting the patient state in coherence with the clinical measurements acquired before stage 2 surgery.

Clinical cases

Patients A, B and C were recruited at the University of Michigan, Ann Arbor, MI, USA; patients D and E at Medical University of South Carolina, Charleston, SC, USA; and patient F recruited at Great Ormond Street Hospital, London, UK. The study was approved by the institutional review board at each institution and consent was obtained from the subject's legal guardian. Patient A had a diagnosis of pulmonary atresia with intact ventricular septum. The stage 1 surgery comprised a central shunt, from the ascending aorta (AoA) to the main pulmonary artery (MPA). This patient was 6 months old at time of cardiovascular magnetic resonance imaging (CMR), with body surface area (BSA) of 0.34 m 2 . Patients B, C, D and F had hypoplastic left heart syndrome. Their stage 1 operations were all 3.5 mm right modified Blalock-Taussig shunts (rmBT). These patients were 4, 5, 5 and 4 months old at time of CMR, with BSA of 0.28, 0.34, 0.30 and 0.27 m 2 respectively. Patient E had tricuspid and pulmonary atresia, and underwent stage 1 surgery of a 4 mm rmBT, with additional left pulmonary artery patch to treat a left pulmonary artery stenosis. This patient was 3 months old at time of CMR, with BSA of 0.26 m 2 .

Clinical measurements

For patients A, B and C, CMR was performed immediately prior to the second surgery under the same general anesthetic, and cardiac catheterization was performed respectively 3 days, 1.5 months and 1 week prior to the second surgery. CMR and cardiac catheterization were performed on the same day, 1.5 months prior to the second surgery for patient D; and 1 month prior for patients E, and F. Flows were quantified by cardiac-gated CMR, and thus they give in the best case cycle-averaged time-varying information. For patients A and D, flows were measured in the branch pulmonary arteries; and for patients B, C, E and F, in the pulmonary veins, for which only time-averages are equal to the corresponding arterial flows. The total pulmonary flow (inflow) was obtained by adding the left and right sided pulmonary flows, and was imposed at the inlet surface (shunt). The mean values are reported in Table 1. Catheterizations were performed in the common atrium and on the venous side, by a pulmonary venous wedge pressure as a surrogate of the pulmonary artery pressure. Pressure differences between left and right sides were either measured or considered negligible (0-1mmHg), except for patient D. The values defined as target pulmonary arterial pressures are reported in Table 1.

The main challenge of this study is to determine outflow boundary conditions for each patientspecific model that are consistent with the clinically measured data: 1) time-average flow distribution between left and right sides (flow split) denoted by fs and 2) pulmonary arterial pressure denoted by P t , which is the target pressure to match, and which location and definition vary for each patient.

Numerical methods

The pre-operative anatomy of each patient was reconstructed from three-dimensional contrast enhanced CMR data using commercial software (Mimics, Materialise, NV, Leuven, Belgium). The vascular region of interest (the distal shunt and the pulmonary arteries extending to their first branches) was selected to include enough branches to avoid numerical instabilities or to largely include the future virtual stage 2 operative anastomosis [START_REF] Corsini | An integrated approach to patient-specific predictive modeling for single ventricle heart palliation. Comput[END_REF][START_REF] Kung | Others: Predictive modeling of the virtual Hemi-Fontan operation for second stage single ventricle palliation: two patient-specific cases[END_REF], without going beyond the branches for which the image resolution would affect the simulation accuracy. A three-dimensional surface was thus obtained, following operations of segmentation and region-growing described in details in [START_REF] Armillotta | Use of rapid prototyping models in the planning of percutaneous pulmonary valved stent implantation[END_REF][START_REF] Schievano | Percutaneous Pulmonary Valve Implantation Based on Rapid Prototyping of Right Ventricular Outflow Tract and Pulmonary Trunk from MR Data[END_REF].

For each of these 3D models, finite element meshes were generated with the commercial software MeshSim (Simmetrix Inc., Clifton Park, MY). Resulting meshes were anisotropically adapted to the flow in several steps [START_REF] Muller | Anisotropic adaptive finite element method for modelling blood flow[END_REF][START_REF] Sahni | Efficient anisotropic adaptive discretization of the cardiovascular system[END_REF]. In order to understand hemodynamics in these single ventricle geometries and reflect the effect of the downstream vascular trees, rigid-wall Navier-Stokes equations were solved in three-dimensional models, coupled to zero-dimensional models (Figure 1). This multidomain approach and its monolithic numerical implementation are described in details in [START_REF] Vignon-Clementel | Outflow boundary conditions for three-dimensional finite element modeling of blood flow and pressure in arteries[END_REF][START_REF] Vignon-Clementel | Outflow boundary conditions for 3D simulations of non-periodic blood flow and pressure fields in deformable arteries[END_REF]. All simulations were performed with a custom version of finite element Navier-Stokes solver, SimVascular (http://simtk.org), assuming blood as an incompressible Newtonian fluid with a density of 1060 kg m -3 and a dynamic viscosity of 0.004 Pa s.

Tuning of reduced parameters

The automatic tuning of outflow boundary condition was introduced in [START_REF] Troianowski | Others: Three-dimensional simulations in Glenn patients: clinically based boundary conditions, hemodynamic results and sensitivity to input data[END_REF] for patient-specific simulations of the stage 2 (Glenn) procedure. This method was developed to obtain 3-element Windkessel parameters containing proximal resistance (R p ) to describe viscous effects into large arteries; capacitance (C) representing the ability of large arteries to store and return blood flow during a cardiac cycle; and distal resistance (R d ) representing viscous effects into small arteries. A first part consisted of iteratively tuning a total resistance at each outlet by coupling 3D Navier-Stokes equations to a resistance with a fixed point method until recovering: 1) the measured transpulmonary gradient between the superior vena cava and the left atrium, and 2) the pulmonary flow split. Here this method is extended in order to take into account the different locations of measured pulmonary arterial pressure: instead of only handling situations where the measured pressure coincides with the pressure at one boundary of the 3D domain (the inlet), the former can be more general as will be explained below. This changes the way the outlet parameters are updated at each tuning iteration, and possibly the algorithm convergence. The framework presented below is the same for steady and pulsatile simulations.

The target flow split (fs) is computed from the time-average right pulmonary (Q R ) and left pulmonary (Q L ) CMR flows as follow:

fs = Q R Q L + Q R (1)
The mean flow rate (Q i ) that is the target to reach through an outlet i is assumed to be proportional to the surface area of this outlet (S i ):

Q i ∝ S i (2)
A power law could also be used if relevant [START_REF] Troianowski | Others: Three-dimensional simulations in Glenn patients: clinically based boundary conditions, hemodynamic results and sensitivity to input data[END_REF]. A target average value for the flow rate at outlet i is thus as a function of the inlet flow rate (Q in ) coming from clinical measurements:

Q i = Q in δ ir fs + (1 -δ ir )(1 -fs) S i j S j δ ij (3)
where δ ir is the kronecker symbol equal to 1 if outlet i is on the right side and 0 if not, and δ ij is the kronecker symbol equal to 1 if outlet j is on the same side as outlet i and 0 if not.

A first value of R i at each outlet is estimated to initialize the algorithm. This serves as the boundary condition for a first 3D Navier-Stokes simulation that is run during enough time steps in order to reach stable state results, typically a few 100 time steps for steady simulations and four cardiac cycles for pulsatile simulations. Note here that even in case of steady inlet flow (so called "steady simulations"), the resulting 3D flow is unsteady due to its complex interaction with the patientspecific geometry. Therefore, we use a transient formulation even for steady boundary conditions.

By post-processing the results over the last stable period, time-averaged (and mean in space) pressure (P 3D i ) and flow rate (Q 3D i ) are computed at each outlet. A flow split (fs 3D , based on the sum of the flows Q 3D i on each side), and a controlled pressure (P 3D c ) that varies according to the measurement location and clinician input, are then calculated. Note that the exact pressure measurement location was unknown, hence this controlled pressure could not be related to the output pressure at a specific location in the 3D domain. Rather, if catheterization was performed on the left pulmonary side, P 3D c is the average pressure over all the outlets of the left side

(P 3D i , i ∈ [0, N ]
for N the number of left pulmonary outlets). Sometimes, there were several measurements done on one or both sides; P 3D c is then defined with the clinical expert as being the maximum or the minimum value over all the corresponding branches of number N , e.g.

P 3D c = min i P 3D i , i ∈ [0, N ] (4)
These different options define the measurement for tuning and control method reported in Table 1. The results are then compared to the target pressure and flow values,

ε P = |P t -P 3D c | P t ε Q = 1 N N i=1 Q 3D i -Q i Q i - 1 N N i=1 Q 3D i -Q i Q i 2 1 2 (5)
where N is the total number of outlets. If convergence is not reached within a certain tolerance, the resistance at each outlet is updated as follows:

R i = P t -(P 3D c -P 3D i ) Q i (6)
Another 3D simulation is run coupled to the new set of resistances R i , where i ∈ [0, N ] for N outlets. This process is repeated until convergence.

When pulsatile tuning simulations are performed, which can be important to better capture nonlinear effects across the 3D domain, resistances are also applied as boundary conditions, while pulsatile flow is imposed at the inlet. The controlled variables (P 3D c , Q 3D i ) are then the timeaveraged values over the last cardiac cycle. The steady tuning simulations were performed for all the six patients, whereas the pulsatile ones only for patients A and E due to time-varying arterial inflow issues. An additional step could be naturally included between the update of equation 6 and the next 3D iterative run, to update and impose more complex downstream reduced models, such as the one presented in the next section. This was not done here because the time-varying target pressure measurements were not available. Hence for each outlet, the 5-parameter reduced model was derived from the total resistance, with a combination of modeling and literature data. This is the object of the next section.

Reduced model parameters for pulsatile simulation

From the pulmonary arteries to the left atrium the vascular tree is carrying blood though the lung where gas exchanges occur. We divide this vasculature in three parts: the arteries downstream of the 3D geometry, the smaller arteries and capillaries, and the veins; representing the viscous effects and the ability of the larger vessels (both arteries and veins) to store and return blood during the cardiac cycle. The resulting parameters described on Figure 1 ) and one capacitance (C ART ) is generated from the Womersley-based impedance of this sub-tree [START_REF] Spilker | Morphometry-Based Impedance Boundary Conditions for Patient-Specific Modeling of Blood Flow in Pulmonary Arteries[END_REF][START_REF] Troianowski | Others: Three-dimensional simulations in Glenn patients: clinically based boundary conditions, hemodynamic results and sensitivity to input data[END_REF]. Hence the proximal to distal pulmonary artery resistance ratio is different for each outlet. The capacitance in the venous tree is 5.7 bigger than in the arterial side. Combining all these relations determines a unique set of reduced parameters for each outlet i based on its total resistance. In a next step, the three dimensional Navier-Stokes equations can be coupled to these reduced models at the outlets for more complete pulsatile simulations. A preoperative overview of the hemodynamics in the pulmonary arteries of each patient is thus generated in order to investigate a palliative surgery [START_REF] Corsini | An integrated approach to patient-specific predictive modeling for single ventricle heart palliation. Comput[END_REF][START_REF] Kung | Others: Predictive modeling of the virtual Hemi-Fontan operation for second stage single ventricle palliation: two patient-specific cases[END_REF].

R i = R ART p + R ART d + R VEIN d (7)

Results

The Results section presents the demonstration of this framework on six stage 1 single ventricle patients, defining their specific simulations parameters and reporting their 3D pressure, flow and wall shear stress. All 3D geometrical models included an inlet which is the distal shunt and a number of pulmonary branch outlets as described in the methods. The inlet surface areas of patients A, B, C, D, E, and F are respectively 3.39 10 -2 , 1.27 10 -1 , 1.83 10 -1 , 6.21 10 -2 , 6.29 10 -2 and 5.93 10 -2 cm 2 . The 3D meshes of patients A, B, C, D, E and F contain respectively 715K, 970K, 1050K, 1500K, 710K and 1540K tetrahedra. Recall that the clinical data are the flow split between the right and left lungs and a pulmonary pressure. For patients A, B, C, D, E, and F, we should recover a flow split equal to 0.64, 0.46, 0.55, 0.46, 0.67 and 0.52 respectively, and a pulmonary arterial pressure equal to 12.8, 12.7, 13.5, 12, 14 and 11 mmHg respectively, depending on the location where the clinical measurements were performed. Table 1 reports the pressure measurement location chosen for the tuning and the controlled value (average, minimum or maximum over branches of the measurement location) that was set by the clinicians. The pressure difference between the right and left pulmonary artery sides was negligible, except for patient D where it was 3 mmHg. Table 1 also reports the tuning results. Steady tuning was run for every patient: the results are very close to the targeted clinical data. For the first patient, pulsatile tuning was also performed for patient A and led to results of 0.64 for the flow split, 12.4 mmHg for the pulmonary pressure and 0.4 mmHg for the pressure difference, all very close to both steady tuning and targeted clinical values. For patient E, who was the only case with significant pressure difference between right and left lungs, results are 0.67, 17 mmHg and 3 mmHg respectively, also very close to both steady tuning and targeted clinical values. As a consequence, the next paragraphs only report the 3D results of the steady tuning. and pressure difference between left and right sides ∆P 3D L/R ). Numbers are for the steady tuning and in parenthesis for the pulsatile tuning.

Pressure distribution

We present in Figure 3 the pressure map for each patient at the end of the tuning. We observe a high peak of pressure on the wall. The maximal pressure for patients A, B, C, D, E and F are respectively 52, 26, 25, 45, 62 and 74 mmHg. The patients A, B and F have homogeneous pressure in both pulmonary arteries respectively, approximately 12, 13 and 11 mmHg. Moreover, patients D and E have a left pulmonary artery (LPA) stenosis which involves a significant pressure loss of approximately 3 mmHg across the stenosis, while homogeneous pressure is found in the right pulmonary artery (RPA). A last remark concerns patient C who has a kink in the RPA close to the shunt. The kink does not generate a significant pressure loss and pressure is homogeneous in both PAs, around 13.5 mmHg.

Velocity patterns

In Figure 4 blood flow patterns for each patient are presented at the end of the tuning. A very complex flow can be observed in the 3D geometrical domain, especially close to the shunt where blood flow is the highest. Highest velocity is obtained in the center of the shunt and is respectively equal to 470, 250, 265, 319, 377 and 570 cm/s. For patients A, B, and F the blood is swirling in the pulmonary arteries close to the shunt and goes towards the outlet smoothly with less complexity. This mixing flow behavior occurs into the main pulmonary artery stump for patients B, C, and D where the anastomosis of the ligated main pulmonary artery is prominent. For patient E, swirling reflects the tortuous geometry on the left side, where there are a constricted anastomosis of the ligated main pulmonary artery, followed by a stenosis.

Wall shear stress maps

In Figure 5 the wall shear stress surface maps are represented for each patient. For all patients the maximal wall shear stress is reached either in the shunt, or where there is a peak of pressure, respectively equal to 73.8, 55.3, 52.9, 66.3, 67.1, and 69.7 Pa.

Moreover, for patients D and E, the stenosis involves a large wall shear stress respectively equal to 54.3 and 61.4 Pa. Regarding patient A, where the main pulmonary artery was ligated the mean wall shear stress is 57.8 Pa, which is high compared to the rest of the domain. In contrast, in the corresponding location of patient E, the mean wall shear stress is close to 0.

Comparing patients B and C, which have a similar shape, the obtained mean wall shear stress is respectively equal to 17 and 40 Pa.

Discussion

Significance of the Results

The proposed methodology to tune reduced model parameters accurately replicates clinical measurements 1, at least within acceptable clinical measurement tolerance. Moreover, we highlight the high complexity of the blood flow patterns in pulmonary arteries in the context of the stage 1 single ventricle procedure. Patients have a BSA between 0.26 and 0.34 m 2 . The variability of the pulmonary flow rate -between 7.5 and 20 cm 3 /s -leads to very high Reynolds numbers, between 1750 and 4150 when computed at the shunt inlet with the maximum velocity.

This leads to values of wall shear stress that are very high. They are in the same order of magnitude as in the aorta of a typical Norwood patient [START_REF] Qian | Computational Hemodynamic Analysis in Congenital Heart Disease: Simulation of the Norwood Procedure[END_REF] (keeping in mind that the shunt is a direct connection from the aorta or brachial artery). In the Glenn circulation [START_REF] Troianowski | Others: Three-dimensional simulations in Glenn patients: clinically based boundary conditions, hemodynamic results and sensitivity to input data[END_REF], a typical wall shear stress is more than 10 times lower. For patients B and C, the main pulmonary artery was ligated, forming a sphere bulging out of the anastomosis, which might foster recirculation and thus low wall shear stress. Yet, the computed wall shear stresses were quite different. The relative differences in wall shear stress between patients may have some bearing on the development of the PAs, however, the stage 1 circulation is only present for a few months before second stage palliative surgery is required. Therefore, these effects will be of a fairly short duration.

The importance of the severe stenosis for the pressure drop and the swirling behavior of the blood flow was also shown. This complexity is induced by the high Reynolds number is in the shunt. For these stenosis cases and patient C with its kinks, computational fluid dynamics is a a complement to clinical measurements. One potential clinical application for this technique is to aid clinicians understanding of the significance of anatomical abnormalities in these complex patients. For example, patient C appears to have a significant geometric 'kink' or restriction in their RPA just distal to the shunt. Only a distal LPA pressure measurement was obtained clinically, so the hemodynamic effect of this lesion was unknown. However, simulation results clearly indicate that the pressure difference between the shunt anastomosis and the PAs are equal on the left and right sides. This was of interest to the clinicians, as the consensus was that judging by the anatomy alone, they would expect a pressure loss across the RPA. Indeed, at stage 2 surgery, patients A, C, D, and E all underwent patch augmentation of their central PAs or proximal LPAs judged on the appearance of the geometry alone. This work suggests that in some cases the clinician's perception of what constitutes a geometric abnormality may actually not result in a hemodynamically significant pressure difference, as in case C. Although, the tortuous course of the PA appears to constitute a stenosis, it is partly an illusion due to the distal PA having a relatively large cross-sectional area.

For patient D, a first attempt was to match the measured flow split of 0.46 and a pressure at both LPA and RPA side equal to 12 mmHg. Unfortunately, targeting the flow split to be 0.46 induced a pressure difference between left and right sides around 6 mmHg which was not acceptable from a clinical point of view. Discussions with clinicians led us to disregard the clinically measured flow split. By imposing the same pressure on average on each side, a flow split of 0.54 (more flow to the right side) was then found, which was clinically acceptable. The differences in flow split would correspond to an error of 0.78 ml/s (i.e. 8 %) in estimating flow with CMR. This example highlights the impact of uncertainty of measurements on the numerical simulations. It could be interesting to investigate more precisely the effect of uncertainties of pressure and flow split measurements on the pressure loss through the stenosis. This patient-specific example was thus challenging and required a close collaboration with clinicians. On the other hand, performing simulations underlined the incoherence between the original clinical measurements, and helped to assess together with the clinical experts which data made more sense to trust.

Methodology Discussion

The tuning methodology introduced in this work consists of coupling Navier-Stokes equations to reduced models (here resistances) at the outlets and running steady or pulsatile simulations. The a priori interest of running pulsatile simulation is the integration of all non linearities due to the 3D geometrical model and obtain a more accurate set of reduced model parameters. However, to achieve periodic stability, pulsatile simulations need to be run over 4 cardiac cycles (around 2000 time steps) and the iterative process is converging in 5 iterations, thus around 10000 time steps. In comparison, 100 time steps for a steady stimulation were enough and the iterative process converges in 4 iterations, thus around 400 time steps were necessary in total. Furthermore, pulsatile tuning for patients A and E were performed and both matched targeted clinical values and were very close to the steady tuning results. Moreover, the reduced model parameters were very close to those obtained by steady tuning. Indeed, for patients A and E the difference between steady and pulsatile tuning are respectively 6.37% and 7.82% in resistances, computed by the formula [START_REF] Bove | Computational fluid dynamics in the evaluation of hemodynamic performance of cavopulmonary connections after the Norwood procedure for hypoplastic left heart syndrome[END_REF]. Note that the same framework could be used to tune total capacitance values for example or even more parameters, if pulsatile pressure measurements were available in synchrony with the flow measurements [START_REF] Pant | A multiscale filtering-based parameter estimation method for patient-specific coarctation simulations in rest and exercise[END_REF].

For sake of brevity, the 5 parameters of each branch for the six patients considered in this study are not presented in the results section. Such reduced models of the pulmonary circulation for each outlet have been incorporated into preoperative closed-loop simulations representative of the patient state before surgery. They have been then used to investigate different surgical options as in [START_REF] Corsini | An integrated approach to patient-specific predictive modeling for single ventricle heart palliation. Comput[END_REF][START_REF] Kung | Others: Predictive modeling of the virtual Hemi-Fontan operation for second stage single ventricle palliation: two patient-specific cases[END_REF], keeping all these outlet blocks the same while changing the 3D part and its connecting inflow. Deriving this 5-parameter model of the arterial and venous pulmonary circulation has lead to more successful matching clinical time-varying tracings in the closed-loop simulations than the previous 3-component Windkessel model that neglected in particular the venous compliance.

A first source of uncertainty comes from the geometrical reconstruction for these young babies. The typical pixel size from 3D CMR angiography is 1 * 1mm 2 . In addition 3D CMR angiography may suffer from signal loss in areas of turbulent or complex flow. Thus including or excluding a pixel in a small area region, such as a shunt stenosis, may lead to errors in the 3D geometry. This kind of uncertainty should be addressed in future work. Where available, 2D CMR black blood imaging sequences were used to check diameters at specific locations since spatial resolution is higher in this sequence. As a final check for each patient, a visual comparison was made between the 3D reconstruction and the 2D cardiac catheter angiography sequences. Patient D had the appearance of luminal narrowing with an obvious distal stenosis on catheter angiographic sequences. The 2D measurements were used to adjust the 3D reconstruction to give an accurate representation of the geometry in this region.

Branching bifurcations were included evenly on both sides, coherently with image resolution. Although a large number of branches increases computational time, cutting the branches before the first bifurcations on each side often leads to numerical instabilities, as was tested for patient A. Stabilizing the numerical scheme is then warranted [START_REF] Arbia | A new approach for the outflow boundary conditions in three-dimensional hemodynamics[END_REF], [START_REF] Moghadam | A comparison of outlet boundary treatments for prevention of backflow divergence with relevance to blood flow simulations[END_REF]. Note that sometimes, as was the case for patient F, the robust convective stabilization was not enough to avoid numerical divergence in the short model. Besides, as cutting branches or stabilization can change flow features close to the boundary, it is necessary to include enough branches to not affect the flow close to the current anastomosis or the one of the next stage virtual surgery. Furthermore, in this work walls of the 3D geometrical models were assumed rigid because data for the heterogeneous assessment of the wall parameters were not known. Moreover, the shunt is a rigid graft and the anastomosis is also quite rigid because of sutures. If time-varying geometry data were available, the elasticity parameters could also be inferred [START_REF] Bertoglio | Sequential parameter estimation for fluid-structure problems: Application to hemodynamics[END_REF]. However, it has been shown that FSI has little effect at rest at least at later stage of the palliation [START_REF] Bazilevs | Computational fluid-structure interaction: methods and application to a total cavopulmonary connection[END_REF].

Regarding the inlet boundary condition, since the velocity profile was not available, a parabolic profile was imposed at the inlet, which is commonly used for blood flow in arteries. We verified on one case that changing the profile to a flat velocity profile did not significantly change the pressure or the velocity fields.

The clinical hemodynamics measurements have uncertainties and it would be interesting to take it into account in a future work [START_REF] Troianowski | Others: Three-dimensional simulations in Glenn patients: clinically based boundary conditions, hemodynamic results and sensitivity to input data[END_REF]. Indeed it would be useful to investigate their effect on the pressure loss through stenosis and on the set of reduced model parameters. A last comment regarding uncertainties of clinical measurements: in Figure 6 simulated velocity vectors for patient D are shown at the same locations as where flow measurements were done in CMR. We highlight the complexity of the blood flow patterns at these locations, which might explain the difficulty to get accurate clinical time-varying velocity measurements by CMR at these locations. The PAs are relatively small, and are receiving flow from the aorta, or one of its branches, at high pressure. CMR flow measurements are typically less reliable in areas of complex flow because unpredictable phase shifts occur. This often leads to signal loss (voiding) and inaccurate flow measurements [START_REF] Oshinski | Turbulent Fluctuation Velocity: The Most Significant Determinant of Signal Loss in Stenotic Vessels[END_REF]. Typically, the PA flow measurements in this location would be underestimated in stage 1 patients. The pulmonary veins may represent a more stable location to measure the pulmonary flow split, however, the physiology of these patients sometimes leads to the pulmonary veins receiving additional blood flow external to the pulmonary arteries.

It is only through the iterative process of comparing simulation results to clinical measurements that suitable strategies for managing measurement uncertainty are reached.

Conclusion

In this work we develop a method to iteratively tune parameters of reduced models of 3D blood flow simulations in order to reflect the effect of the downstream vasculature, taking into account clinical measurements of stage 1 single ventricle patients before they underwent the next palliative operation. This method was effectively demonstrated on six patients. It can be useful in other multi-branched blood flow simulations for which pressure or flow measurements are taken at various locations, and chosen to be matched as a minimum, maximum, average or other combinations of values in a given region. Simulations in these stage 1 cases highlighted the complexity of the blood flow patterns into the shunt anastomosis and the pulmonary arteries, as well as the very specific pressure map. The wall shear stress was found to be very high in magnitude, but varying in space differently according to the patients anatomy, not necessarily leading to a low value in the ligated main pulmonary stump. Simulations also assessed the effects of anatomical kinks or stenoses, demonstrating the potential use of the technique: clinicians do not always know whether a lesion is hemodynamically significant. These types of simulations could be helpful when planning whether to intervene on a PA anatomical problem. It involves a lot more surgery to intervene on a PA stenosis so if clinicians better understood the difference between a significant stenosis and an anatomical anomaly that does not have an effect, it would save time and danger intervening when it was not necessary. This methodology is moreover helpful for virtual surgical planning testing different surgical approaches [START_REF] Corsini | An integrated approach to patient-specific predictive modeling for single ventricle heart palliation. Comput[END_REF][START_REF] Kung | Others: Predictive modeling of the virtual Hemi-Fontan operation for second stage single ventricle palliation: two patient-specific cases[END_REF]. It is also shown to be a tool for clinicians to assess coherence of clinical measurements : the iterative approach between comparing clinical data and the simulation results emerged as efficient to understand the limitations in clinical measurements and decide upon a strategy to replicate the clinical case that makes sense to clinicians. This process raises the issue of the impact of uncertainty quantification of the clinical measurements.

Finally, this first work on single ventricle shunt anastomoses at stage 1 could be used to compare hemodynamics of different stages of the Fontan palliation.

  need thus to be determined. Assuming that pressures return periodically to the same values in the 5-parameter reduced model, the proximal arterial resistance R ART p , the distal arterial resistance R ART d and the venous resistance R VEIN d can be summed to recover the total resistance at outlet i,

Fig. 1

 1 Fig.1Multiscale set-up, with the 3D domain, its inlet face on which the inflow Q in is prescribed, and the distal pulmonary circulation for each outlet that all merge in the single atrium defined by its pressure Pat. At each outlet (dashed rectangle), the pulmonary arteries, capillaries and veins are represented either by a resistance for tuning, or by a more complete 5-parameter reduced model for final representation.

3. 1

 1 Figure 2 shows the reconstructed 3D geometries of patients A, B, C, D, E, and F.

Fig. 2

 2 Fig. 2 3D geometrical model of patients A, B, C, D, E, F.

Fig. 3

 3 Fig. 3 Pressure maps for patients A-F. Maximum pressure is equal to 25 mmHg on the color scale, even if the real maximal pressure is larger. 1 mmHg = 133.3 Pa

Fig. 4

 4 Fig. 4 Streamlines of patients A-F colored by velocity magnitude.

Fig. 5

 5 Fig. 5 Wall shear stress map of patients A-F in Pa. Maximum color scale is fixed at 50 Pa but the real maximum value is larger.

Fig. 6

 6 Fig. 6 Velocity vectors (coloured by their magnitudes) from the 3D simulation shown at the locations where the CMR flow measurements were performed.

Table 1

 1

	Patient	A	B	C	D	E	F
	Age (months)	6	4	5	5	3	4
	BSA (m 2 )	0.34	0.28	0.34	0.30	0.26	0.27
	Inflow (cm 3 /s)	7.5	12	20	9.7	12	16.6
	Remax	2300	1750	2900	2500	3000	4150
		fs	0.64	0.46	0.55	0.46	0.67	0.52
	PA pressure (mmHg)	12.8	12.7	13.5	12	17	11
	∆P L/R (mmHg)	negl.	negl.	negl.	negl.	3	negl.
	Meas. location for tuning	both	LPA	LPA	both	RPA	both
	Control method	average	average	average	average	min	max
		fs 3D	0.64 (0.64)	0.46	0.55	0.46	0.66 (0.67)	0.52
	P 3D c	(mmHg)	12.3 (12.4)	12.7	13.4	12.0	17.0 (17)	11.0
	∆P 3D L/R (mmHg)	0.4 (0.4)	0.7	0.5	6.0	3.0 (3)	0.2

Clinical measurements, numerical method and results. For each patient, clinical data (age, BSA, average pulmonary inflow, deduced Reynolds number Remax, flow split fs, target pulmonary artery pressure value, negligible or not left/right pulmonary pressure difference), tuning set-up (left or right pulmonary side, and control method to define P 3D c ), and numerical results in red (obtained flow split fs 3D , controlled pressure P 3D c

Acknowledgements This work was supported by the Leducq Foundation as part of the Transatlantic Network of Excellence for Cardiovascular Research "Multi-scale modeling of single ventricle hearts for clinical decision support", and by the Associated team Cardio INRIA grant.

Conflict of interest: None.


Hassan Maatouk

Xavier Bay

A New Rejection Sampling Method for Truncated Multivariate Gaussian Random Variables Restricted to Convex Sets

Statistical researchers have shown increasing interest in generating truncated multivariate normal distributions. In this paper, we only assume that the acceptance region is convex and we focus on rejection sampling. We propose a new algorithm that outperforms crude rejection method for the simulation of truncated multivariate Gaussian random variables. The proposed algorithm is based on a generalization of Von Neumann's rejection technique which requires the determination of the mode of the truncated multivariate density function. We provide a theoretical upper bound for the ratio of the target probability density function over the proposal probability density function. The simulation results show that the method is especially efficient when the probability of the multivariate normal distribution of being inside the acceptance region is low.

Introduction

The need for simulation of truncated multivariate normal distributions appears in many fields, like Bayesian inference for truncated parameter space [START_REF] Gelfand | Bayesian Analysis of Constrained Parameter and Truncated Data Problems Using Gibbs Sampling[END_REF] and [START_REF] Geweke | Exact Inference in the Inequality Constrained Normal Linear Regression Model[END_REF], Gaussian processes for computer experiments subject to inequality constraints [START_REF] Da Veiga | Gaussian process modeling with inequality constraints[END_REF], [START_REF] Emery | Simulating Large Gaussian Random Vectors Subject to Inequality Constraints by Gibbs Sampling[END_REF] and [START_REF] Freulon | Conditioning a Gaussian model with inequalities[END_REF] and regression models with linear constraints (see e.g. [START_REF] Geweke | Efficient Simulation from the Multivariate Normal and Student-t Distributions Subject to Linear Constraints and the Evaluation of Constraint Probabilities[END_REF] and [START_REF] Yu | Efficient Algorithms for Generating Truncated Multivariate Normal Distributions[END_REF]). In general, we have two types of methods. The first ones are based on Markov chain Monte Carlo (McMC) simulation [START_REF] Casella | Explaining the Gibbs sampler[END_REF], [START_REF] Laud | Sampling Some Truncated Distributions Via Rejection Algorithms[END_REF] and [START_REF] Robert | Monte Carlo Statistical Methods[END_REF], as the Gibbs sampling [START_REF] Breslaw | Random sampling from a truncated multivariate normal distribution[END_REF], [START_REF] Geweke | Efficient Simulation from the Multivariate Normal and Student-t Distributions Subject to Linear Constraints and the Evaluation of Constraint Probabilities[END_REF], [START_REF] Griffiths | A Gibbs sampler for the parameters of a truncated multivariate normal distribution[END_REF], [START_REF] Kotecha | Gibbs sampling approach for generation of truncated multivariate Gaussian random variables[END_REF], [START_REF] Li | Efficient sampling method for truncated multivariate normal and Student t-distribution subject to linear inequality constraints[END_REF], [START_REF] Robert | Simulation of truncated normal variables[END_REF] and [START_REF] Rodriguez-Yam | Efficient Gibbs Sampling of Truncated Multivariate Normal with Application to Constrained Linear Regression[END_REF]. They provide samples from an approximate distribution which converges asymptotically to the true one. The second ones are exact simulation methods based on rejection sampling (Von Neumann [START_REF] Von Neumann | Various Techniques Used in Connection with Random Digits[END_REF]) and its extensions, [START_REF] Devroye | Non-Uniform Random Variate Generation[END_REF], [START_REF] Hörmann | Automatic Nonuniform Random Variate Generation[END_REF] and [START_REF] Laud | Sampling Some Truncated Distributions Via Rejection Algorithms[END_REF]. In this paper, we focus on the second type of methods.

Hassan Maatouk • Xavier Bay Ecole Nationale Supérieure des Mines de St-Etienne, 158 Cours Fauriel, Saint-Etienne, France e-mail: hassan.maatouk@mines-stetienne.fr, e-mail: bay@emse.fr Recently, researchers in statistics have used an adaptive rejection technique with Gibbs sampling [START_REF] Geweke | Efficient Simulation from the Multivariate Normal and Student-t Distributions Subject to Linear Constraints and the Evaluation of Constraint Probabilities[END_REF], [START_REF] Gilks | Adaptive Rejection Sampling for Gibbs Sampling[END_REF], [START_REF] Martino | An adaptive accept/reject sampling algorithm for posterior probability distributions[END_REF], [START_REF] Martino | A novel rejection sampling scheme for posterior probability distributions[END_REF] and [START_REF] Robert | Simulation of truncated normal variables[END_REF]. Let us mention that in one dimension rejection sampling with a high acceptance rate has been developed by Robert [START_REF] Robert | Simulation of truncated normal variables[END_REF], and Geweke [START_REF] Geweke | Efficient Simulation from the Multivariate Normal and Student-t Distributions Subject to Linear Constraints and the Evaluation of Constraint Probabilities[END_REF]. In [START_REF] Robert | Simulation of truncated normal variables[END_REF] Robert developed simulation algorithms for one-sided and two-sided truncated normal distributions. Its rejection algorithm is based on the uniform distribution. The multidimensional case where the acceptance region is a convex subset of R d is based on the same algorithm using the Gibbs sampling to reduce the simulation problem to a sequence of one-dimensional simulations. In this case, the method requires the determination of slices of the convex acceptance region. Also, Geweke [START_REF] Geweke | Efficient Simulation from the Multivariate Normal and Student-t Distributions Subject to Linear Constraints and the Evaluation of Constraint Probabilities[END_REF] proposed an exponential rejection sampling to simulate a truncated normal variable. The multidimensional case is deduced by using the Gibbs algorithm. In one-dimension, Chopin [START_REF] Chopin | Fast Simulation of Truncated Gaussian Distributions[END_REF] designed an algorithm that is computationally faster than alternative algorithms. A multidimensional rejection sampling to simulate a truncated Gaussian vector outside arbitrary ellipsoids has been developed by Ellis and Maitra [START_REF] Ellis | Multivariate Gaussian Simulation Outside Arbitrary Ellipsoids[END_REF]. For higher dimensions, Philippe and Robert [START_REF] Philippe | Perfect simulation of positive Gaussian distributions[END_REF] developed a simulation method of a Gaussian distribution restricted to positive quadrants. Also, Botts [START_REF] Botts | An accept-reject algorithm for the positive multivariate normal distribution[END_REF] improves an accept-reject algorithm to simulate positive multivariate normal distributions.

In this article, we develop a new rejection technique to simulate a truncated multivariate normal distribution restricted to any convex subset of R d . The method only requires the determination of the mode of the probability density function (pdf) restricted to the acceptance region. We provide a theoretical upper bound for the ratio of the target probability density function over the proposal probability density function.

The article is organized as follows. In Section 2, we recall the rejection method. Then, we present our new method, called rejection sampling from the mode (RSM) and we give the main theoretical results and the associated algorithm. In Section 3, we compare RSM with existing rejection algorithms.

Multivariate Normal Distribution

The General Rejection Method

Let f be a probability density function (pdf) defined on R d . Von Neumann [START_REF] Von Neumann | Various Techniques Used in Connection with Random Digits[END_REF] proposed the rejection method, using the notion of dominating density function. Suppose that g is another density function close to f such that for some finite constant c ≥ 1, called rejection constant,

f (x) ≤ cg(x), x ∈ R d . (1)
The acceptance/rejection method is an algorithm for generating random samples from f by drawing from the proposal pdf g and the uniform distribution.

Theorem 1 (Rejection Sampling Algorithm, Von Neumann [START_REF] Von Neumann | Various Techniques Used in Connection with Random Digits[END_REF]). Suppose that f and g are two pdfs such that f (x) ≤ cg(x) for all x in the support of f . Then the random variable X resulting from the following algorithm is distributed according to f . 1. Generate X with density g.

Generate U uniformly on

[0, 1]. If cg(X)U ≤ f (X), accept X; otherwise, go back to step 1.
Furthermore it can be shown that the acceptance rate is equal to 1/c. In practice it is crucial to get a small c.

Notice that the rejection sampling algorithm is immediately extended to pseudodensity functions (i.e. positive function with finite integral), avoiding the computation of normalizing constant.

Corollary 1. Let C be a subset of R d and f and g be two pseudo-density functions on

C such that f (x) ≤ k g(x). Then the algorithm in Lemma 1 is still valid if the inequality condition cg(X)U ≤ f (X) is replaced by k g(X)U ≤ f (X).
(2)

The rejection constant is c = k C g(t)dt C f (t)dt . Proof. We have f (x) ≤ k g(x)
, and so

f (x) = f (x) C f (t)dt ≤ c g(x) C g(t)dt = cg(x), (3)
with c = k C g(t)dt C f (t)dt . The condition cg(X)U ≤ f (X) is equivalent to k g(X)U ≤ f (X). ⊓ ⊔

Rejection Sampling from the Mode

Suppose that X has multivariate normal distribution with probability density function:

f (x | µ, Σ ) = 1 (2π) d/2 | Σ | 1/2 exp - 1 2 (x -µ) ⊤ Σ -1 (x -µ) , x ∈ R d (4)
where µ = E[X] and Σ is the covariance matrix, assumed to be invertible.

We consider a convex subset C of R d representing the acceptance region. We assume that µ does not belongs to C , which is a hard case for crude rejection sampling. Furthermore, as explained in Remark 1 (see below) the proposed method is not different from crude rejection sampling if µ ∈ C . Without loss of generality, let µ = 0. Our aim is to simulate the multivariate normal distribution X restricted to the convex set C . The idea is twofold. Firstly, we determine the mode µ * correspond- ing to the maximum of the probability density function f restricted to C . It is the solution of the following convex optimization problem:

µ * = argmin x∈C 1 2 x ⊤ Σ -1 x. ( 5
)
Secondly, let g be the pdf obtained from f by shifting the center to µ * :

g(x | µ * , Σ ) = 1 (2π) d/2 | Σ | 1/2 exp - 1 2 (x -µ * ) ⊤ Σ -1 (x -µ * ) . (6)
Then we prove in the next two theorems that g can be used as a proposal pdf for rejection sampling on C , and we derive the optimal constant. Theorem 2. Let f and g be the pseudo-density functions defined as

f (x) = f (x | 0, Σ )1 x∈C and g(x) = g(x | µ * , Σ )1 x∈C ,
where f and g are defined in ( 4) and [START_REF] Da Veiga | Gaussian process modeling with inequality constraints[END_REF]. Then there exists k such that f (x) ≤ k g(x) for all x in C and the smallest value of k is:

k * = exp - 1 2 (µ * ) ⊤ Σ -1 µ * . (7)
Proof. Let us start with the one-dimensional case. Without loss of generality, we suppose that C = [µ * , +∞[, where µ * is positive and Σ = σ 2 . In this case, the con- dition f (x) ≤ k g is written

∀x ≥ µ * , e -x 2 2σ 2 ≤ ke - (x-µ * ) 2 2σ 2 , (8)
and so

k * = e (µ * ) 2 2σ 2 max x≥µ * e -xµ * σ 2 = e (µ * ) 2 2σ 2 e -min x≥µ * xµ * σ 2 = e - (µ * ) 2 2σ 2 . (9)
In the multidimensional case, we have

k * = max x∈C e 1 2 (µ * ) ⊤ Σ -1 µ * -x ⊤ Σ -1 µ * . Since µ * ∈ C , we only need to show that ∀x ∈ C , x ⊤ Σ -1 µ * ≥ (µ * ) ⊤ Σ -1 µ * . ( 10
)
The angle between the gradient vector Σ -1 µ * of the function 1 2 x ⊤ Σ -1 x at the mode µ * and the vector (x -µ * ) is acute for all x in C since C is convex (see Figure 1). Therefore, (x -µ * ) ⊤ Σ -1 µ * is non-negative for all x in C .

⊓ ⊔

By now, we can write algorithm 1 as follows:

Theorem 3 (RSM Algorithm). Let f and g be the pseudo-density functions defined as Fig. 1 Scalar product between the gradient vector Σ -1 µ * of the function 1 2 x ⊤ Σ -1 x at µ * and (x -µ * ). The blue lines are the level curves of the function x

→ 1 2 x ⊤ Σ -1 x. f (x) = f (x | 0, Σ )1 x∈C and g(x) = g(x | µ * , Σ )1 x∈C ,
where f and g are defined by ( 4), ( 5) and ( 6). Then the random vector X resulting from the following algorithm is distributed accorded to f .

1. Generate X with pseudo-density g.

Generate U uniformly on

[0, 1]. If U ≤ exp (µ * ) ⊤ Σ -1 µ * -X ⊤ Σ -1 µ * , accept X; otherwise go back to step 1.
Proof. We applied Corollary 1 with the optimal constant k * of Theorem 2. The inequality condition (2) is equivalent to

U ≤ e 1 2 (µ * ) ⊤ Σ -1 µ * e -1 2 X ⊤ Σ -1 X e 1 2 (X-µ * ) ⊤ Σ -1 (X-µ * ) , (11)
which is equivalent to

U ≤ exp (µ * ) ⊤ Σ -1 µ * -X ⊤ Σ -1 µ * . (12)
⊓ ⊔ Remark 1. In practice, we use a crude rejection method to simulate X with pseudodensity g in the RSM algorithm. So if µ ∈ C , RSM degenerates to crude rejection sampling since µ * = µ and f = g. Therefore, the method RSM can be seen as a generalization of rejection sampling.

Remark 2. Our method requires only the maximum likelihood of the pdf restricted to the acceptance region. It is the mode of the truncated multivariate normal distribution. The numerical calculation of it is a standard problem in the minimization of positive quadratic forms subject to linear constraints, see e.g. [START_REF] Boyd | Convex Optimization[END_REF].

Performance Comparisons

To investigate the performance of the RSM algorithm, we consider a zero-mean bivariate Gaussian random vector x with covariance matrix Σ , equal to 4 2.5 2.5 2 .

Assume that the convex set C ∈ R 2 is defined by the inequality constraints:

-10 ≤ x 2 ≤ 0 and x 1 ≥ -15, 5x 1 -x 2 + 15 ≤ 0. ( 13
)
It is the acceptance region used in Figures 2 and3. By minimizing a quadratic form subject to linear constraints, we find the mode

µ * = arg min x∈C 1 2 x ⊤ Σ -1 x ≈ (-3.4, -2.0), (14)
and then we compare crude rejection sampling to RSM. In Figure 2, we use crude rejection sampling in 2000 simulations of a N (0, Σ ). Given the number of points in C (red points), it is clear that the algorithm is not efficient. The reason is that the mean of the bivariate normal distribution is outside the acceptance region. In Figure 3, we first simulate from the shifted distribution centered at the mode with same covariance matrix Σ (step one of the RSM algorithm). Now in the second step of the RSM algorithm, we have two types of points (red and black ones) in the convex set C . The black points are in C but do not respect the inequality constraint [START_REF] Geweke | Exact Inference in the Inequality Constrained Normal Linear Regression Model[END_REF]. The red points are in C , and respect [START_REF] Geweke | Exact Inference in the Inequality Constrained Normal Linear Regression Model[END_REF]. We observe that RSM outperforms crude rejection sampling, with acceptance rate of 21% against 3%. The performance of the method appears when the probability to be inside the acceptance region is low. In Table 1, we consider the one dimensional case d = 1 and we only change the position of µ -. From the last column, we observe that our algorithm outperforms crude rejection sampling. For instance, the proposed algorithm is approximately 14000 times faster than the crude rejection sampling when the acceptance region is [4.5, +∞[. Note also that the acceptance rate remains stable for large µ -(near 10%) for the RSM method whereas it decreases rapidly to zero for crude rejection sampling.

Table 2 Comparison of average acceptance rate between Robert's method [START_REF] Robert | Simulation of truncated normal variables[END_REF] and RSM under the variability of the distance between µ -and µ + . The acceptance region is C = [µ -, µ + ], where µ -is fixed to 1. Robert [START_REF] Robert | Simulation of truncated normal variables[END_REF] also proposed a rejection sampling method in the one dimensional case. To compare the acceptance rates of RSM with Robert's method, we consider a standard normal variable truncated between µ -and µ + with µ -fixed to 1. In Robert's method, the average acceptance rate is high when the acceptance interval is small (see Table 2.2 in [START_REF] Robert | Simulation of truncated normal variables[END_REF]). In the proposed algorithm, simulating from shifted distributions (first step in the RSM algorithm) leads to the fact that the average acceptance rate is more important when the acceptance interval is large. As expected, the performance of the proposed algorithm appears when we have a large gap between µ -and µ + , as shown in Table 2. Now we investigate the influence of the problem dimension d. We simulate a standard multivariate normal distribution X restricted to C = [µ -, +∞[ d , where µ - is chosen such that P(X ∈ C ) = 0.01. The mean of the multivariate normal distribution is outside the acceptance region. Simulation of truncated normal distributions in multidimensional cases is a difficult problem for rejection algorithms. From Table 3, we can remark that when the dimension increases, the parameter µ -tends to zero. Hence, the mode µ * = (µ -, • • • , µ -) tends to the zero-mean of the Gaussian vector X. And so, the acceptance rate of the proposed method converges to the acceptance rate of the crude rejection sampling. As an additional example when µ -is fixed to 1.35 and d = 5, the RSM algorithm is 135 times faster than the crude rejection sampling. In that case, the probability of the normal distribution X being inside the acceptance region is low.

µ + -µ -Robert

Conclusion

In this paper, we develop a new rejection technique, called RSM, to simulate a truncated multivariate normal distribution restricted to any convex set. Our method only requires to find the mode of the target probability density function restricted to the convex acceptance region. The proposal density function in the RSM algorithm is the shifted target distribution centered at the mode. We provide a theoretical formula of the optimal constant such that the proposal density function is as close as possible to the target density. Note that the RSM algorithm is easy to implement. An illustrative example to compare RSM with crude rejection sampling is included. The simulation results show that using rejection sampling from the mode is more efficient than crude rejection sampling. Comparisons with Robert's method in the one dimensional case is discussed. The RSM method outperforms Robert's method when the acceptance interval is large and the probability of the normal distribution to be inside is low. The proposed rejection method has been applied in the case where the acceptance region is a convex subset of R d , and can be extended to non-convex regions by using the convex hull.

Fig. 2

 2 Fig. 2 Crude rejection sampling using 2000 simulations. The acceptance rate is 3%.

Fig. 3

 3 Fig. 3 Rejection sampling from the mode using 2000 simulations. The acceptance rate is 21%.

Table 1

 1 Comparison between crude rejection sampling and RSM when the probability to be inside the acceptance region becomes low. The acceptance region is C = [µ -, +∞[.

	µ -Average of acceptance rate with Average of acceptance rate with Gain
		crude rejection sampling (%)	RSM (%)
	0.5	30.8	34.9	1.1
	1	15.8	26.2	1.6
	1.5	6.7	20.5	3.0
	2	2.2	16.8	7.4
	2.5	0.6	14.2	23.1
	3	0.1	12.2	92.0
	3.5	0.0	10.6	455.6
	4	0.0	9.3	2936.7
	4.5	0.0	8.4	14166.0

Table 3

 3 Comparison between crude rejection sampling and RSM with respect to the dimension d. The acceptance region is C = [µ -, +∞[ d .

	Dimension µ -Average of acceptance rate with Average of acceptance rate with Gain
	d		crude rejection sampling (%)	RSM (%)
	1	2.33	1.0	15.0	15.0
	2	1.29	1.0	5.2	5.2
	3	0.79	1.0	2.5	2.5
	4	0.48	1.0	1.5	1.5
	5	0.25	1.0	1.2	1.2

Acknowledgements

The authors wish to thank Olivier Roustant, Laurence Grammont and Yann Richet for helpful discussions, as well as the participants of MCQMC conference.


P Baracchini
email: paul.baracchini@u-bordeaux.fr

C Guillebaud

F X Kromm

A Catapano

M Montemurro

H Wargnier

ARCHITECTURE AND MATERIALS SELECTION IN MULTI-MATERIALS DESIGN

Keywords: Multi-materials, Architectured materials, Materials selection, Optimisation

The design process involving both the architecture and the materials represents an hard task mainly due to the high number of potential configurations, thus requiring firstly the development of new and more rigorous approaches but also the development of new tools. To this purpose, we present in this work a new strategy for the design of architectured materials. Such a strategy relies on one hand on the construction of some databases for the selection of both geometrical patterns and materials, and on the other hand on the use of well-known analytical models to describe the physical behaviour of the multi-material. In order to prove its effectiveness, we apply our strategy to the problem of the least-weight design of a multilayer plate that has to meet thermal, electrical and mechanical requirements. Moreover, we use a genetic algorithm, as a numerical tool, to perform the solution search for our problem. Numerical results show that we can obtain optimum configurations characterised by a weight saving up to 59% keeping the same (or even superior) thermal, electrical and stiffness properties than those of a monolithic reference plate.

Introduction

Designing an industrial product consists in determining all its characteristics: number of parts, geometry, dimensions, constitutive materials, manufacturing process, joining process, etc. Several studies were conducted in order to elaborate a strategy putting all these characteristics into the same design process by developing, for instance, different analytical methods like functional analysis, concurrent engineering or TRIZ. Among these contributions, materials selection methods are based on three different approaches [START_REF] Ashby | Selection strategies for materials and processes[END_REF][START_REF] Ashby | Materials selection in mechanical design[END_REF]:

(1) the free search method which explores the whole set of solutions;

(2) the questionnaire strategy which guides the designer by asking appropriate questions;

(3) the knowledge based system which uses the results of previous experience.

The fundamental principle of materials selection based on free searching lies on the quantitative evaluation of the performance of a material regarding the functions of the product. As the best solution is picked from a database, it offers better possibilities for innovation than the others. However, the use of performance index isn't so simple when designing a multi-material.

In accordance with the definitions proposed by Ashby, Bréchet or Kromm [START_REF] Ashby | Designing hybrid materials[END_REF][START_REF] Brechet | Architectured materials: expanding materials space[END_REF][START_REF] Kromm | Definition of a multimaterials design method[END_REF] a multimaterial or an architectured material is considered to be the association of one or several materials disposed according to a predefined architecture such that a representative volume element has at least one dimension that is very small compared with the dimensions of the entire structure.

The parameters that the designer must define in the design of a multi-material are:

the components, that are usually materials, or be semi-products (this is the case, for example, of multilayer structures or stratified composites); the volume fractions of the components; the architecture and morphology of the components, i.e. their spatial disposition; coupling modes between the components, especially the nature of the interfaces and their behaviour. Previous studies have been carried out to find a way to determine some of these parameters. For example, a focus on the material requirements allows the division of the constraints in two sets of requirements, separating the incompatible ones [START_REF] Wargnier | Study of the compatibility between criteria in a set of materials requirements. Application to a machine tool frame[END_REF]. Moreover, when the architecture of the multi-material has been chosen, the problem is reduced to the search of a proper combination of materials, but even in this case the number of configurations remains very high. To this purpose, in [START_REF] Giaccobi | Filtration in materials selection and multi-materials design[END_REF] the authors proposed a filtration method to determine all possible combinations of materials satisfying the constraints imposed to the problem. In this way the number of potential solutions, i.e. those configurations that are candidate to be solutions for the problem at hand is drastically reduced, thus, simplifying the solution search process. The aim of this paper is to show how the architecture of the multi-material can be selected by the designer. In a first paragraph, the general principle is described. Then, as this method is based on a free searching approach, an architecture database has been built. Finally, the numerical method for the choice of the optimal architecture and components is explained and validated on an example.

Basis for multi-material architecture selection

General principle of the selection method

The selection of the best material for an application consists in determining the optimal value of a combination of material properties called performance index. This indicator is derived from the expression of the performance of the product (mass or cost for example) in which all the parameters except the materials properties (i.e. geometrical or functional parameters) are fixed. There are several ways to define the performance of architecture. For example, the improvement on a property (or combination of properties), i.e. the gain, could be calculated comparing it with a law of mixture with the same components and volume fraction. This principle is close to the shape factor [START_REF] Weaver | The optimal selection of material and section-shape[END_REF][START_REF] Pasini | Shape transformers for materials and shape selection of lightweight beams[END_REF] that is used to select the shape of the cross section of a beam. Another method consists in using shape optimisation algorithms like level set for example [START_REF] Laszczyk | Homogenization and topological optimization of architectured panels[END_REF]. With this method, the architecture isn't picked from a database, but it evolves at each step of the calculation. The drawbacks of this method are that the final morphology depends greatly on the initial one, and that the computational costs times are very expensive. The method developed in this study consists in selecting the architecture from a database, as shown in figure 1. As the homogenised properties of this architectured material depend on the The selection in a database allows a quick comparison of very different morphologies. Then, after a first ranking of the candidates, the parameters of the chosen architecture can be optimised thanks to classical topology optimisation methods in a final step. Before starting the selection process, the architecture database has to be built, taking into account the important notion of hierarchy.

Architecture database

The database that must be created has to be representative of the various possibilities of the multi-materials, but doesn't have to be exhaustive because it just aims at illustrating the selection method and mustn't lead to important calculation times. Trying to make a collection of the most classical architectured materials, it appears fundamental to create a hierarchy in the database. Indeed, multi-material morphologies can sometimes appear quite complex, but they can be considered as a combination of elementary patterns [START_REF] Dunlop | Multilevel architectures in natural materials[END_REF][START_REF] Bollen | Multifunctional architecture materials for electromagnetic absorption[END_REF], so, the database has to be organized as a function of the length scale. In a first time, at macroscopic scale, it is supposed here that all the multi-materials are multilayered, so the only morphology that is considered in this category is a division of the material through the thickness. Then, each layer can be filled with a predefined morphology. The possible elementary patterns have been separated in different types as illustrated in figure 2: 1) monolithic material;

2) composite material: matrix reinforced by particles, short or continuous fibers; 3) cellular material like foams or honeycombs; 4) functional patterns, i.e. segmented morphologies that allow a fluid flow circulation.

Homogenisation models

After the presentation of the most representative patterns, the corresponding homogenisation models have to be identified. The analytical homogenisation model efficiency is at the crossroads between the right choice of the number of geometrical parameters considered and their usability. A classification of the homogenisation models is for example made for the composite materials [START_REF] Mottram | Thermal transport properties[END_REF] following the parameters that are taken into account (morphology of reinforcement, statistical repartition in the matrix…). The chosen models for this study are The interfaces properties are not taken into account as a design parameter, the layers are supposed to be perfectly bonded.

Numerical method for architecture and components selection

Definition of a candidate

In this study, 22 parameters have to be determined for each layer, with 8 to 100 levels depending of parameter. The available range of value associated to each parameter, and the step that is used to discretise these variables is an important way of the design. Indeed, the interval of variation must on one hand take into account the physical or technological limitations. But, the chosen variation step must on the other hand be high enough to limit the number of candidates, but low enough to avoid the elimination of solutions because of a rough dimensioning.

Evaluation and ranking of the solutions

The solution space is too big to allow a systematic screening of the solution space because of the calculation time into play. The numerical method for the selection of the architecture and components has to deal with a lot of different parameters to evaluate the homogenised properties and performances. Some of these parameters are quantitative (volume fractions, particle dimensions…), but others are qualitative (type of architecture, constitutive materials…) so different kinds of variable have to be manipulated. In this case, previous studies showed that genetic algorithms can give optimum results in reduced times [START_REF] Leite | Optimal design of architectured sandwich panels for multifunctionnal properties[END_REF][START_REF] Chakraborti | Tailor-made material design:An evolutinary approach using multi-objective genetic algorithms[END_REF], so we chosed this numerical method.

Case study

Description of the design problem

The developed method has to be, firstly, validated on a basic case study. This example concerns the weight minimisation of a rectangular bi-layer plate, satisfying some design constraints on the in-plane tensile stiffness value along one direction (case 1) or two orthogonal directions (case 2), on the maximum temperature on the bottom face. This plate is submitted to a given transverse heat flow and an imposed in-plane electrical resistance, as illustrated in figure 3. The transverse heat flow is imposed on the lower side of the plate with a value equal to 900 W/m 2 . The thermal modelisation takes into account the heat conduction through the plate, and the convection on the upper face. As the convection coefficient depends not only on the air flow but on the geometry of the surface, it is attributed a value of 18 W/Km 2 for a plane surface, but has to be calculated when the upper layer of the material is constituted of fins [START_REF] Khan | Modeling of cylindrical pin-fin heat sinks for electronic packaging[END_REF][START_REF] Teertstra | Analytical forced convection modeling of plate fin heat sinks[END_REF]. Finally, electrical resistance along the length direction of the plate is required smaller than 2 mΩ.

Optimal solutions for the problem

To solve the problem, we used the genetic algorithm BIANCA [START_REF] Montemurro | A two-level procedure for the global optimum design of composite modular structures -application to the design of an aircraft wing. Part 2: numerical aspects and examples[END_REF][START_REF] Montemurro | Optimal design of advanced engineering modular systems through a new genetic approach[END_REF] which can handle, within the same problem, variables of different nature: continuous, discrete, scattered and abstract. Depending on the various mechanical loads imposed, i.e. case 1 or 2, the following results were obtained. For comparison of the results, reference is made to a 2 mm thick plate made of aluminium alloy dimensions 250 mm x 100 mm. From the results of table 1, it can be noticed that the genetic algorithm shows good matching in the architecture choice responsive to a mechanical load change. Moreover, these architecture/materials couple selection were obtained in a reasonable time (a few minutes).

Temperature evaluation

Heat flow direction (conduction + convection) Imposed stiffness

Conclusions

A method allowing the selection of the architecture and components of a multi-material has been developed. Thanks to a database presenting the most classical architectures and to a material database, the genetic algorithm used in this work has been able to propose optimised solutions in a very short time with respect to the number of potential solutions. This method has been validated on a schematic example lying on mechanical, thermal and electrical constraints. The results illustrate the adaptability of the multi-materials, giving different optimal solutions when constraints were changed. In order to make this selection more precise, it would be interesting to establish filtration criteria for architectures and components, so that the solution space would be narrowed. With an efficient method, the reduced number of candidates would allow a complete screening with smaller variation steps for more precise selection.

Layer

Figure 1 .

 1 Figure 1. General method for architecture selection

Figure 2 .

 2 Figure 2. Hierarchical database of elementary patterns

Figure 3 .

 3 Definition of the design constraints for (a) mechanical and (b) thermal functions

Table 1 .

 1 Results of the optimisation calculations

	top 1	Layer bot 2	% Weight gain