Cardiovascular Modeling With Adapted Parametric Inference

Computational modeling of the cardiovascular system, promoted by the advance of fluid-structure interaction numerical methods, has made great progress towards the development of patient-specific numerical aids to diagnosis, risk prediction, intervention and clinical treatment. Nevertheless, the reliability of these models is inevitably impacted by rough modeling assumptions. A strong in-tegration of patient-specific data into numerical modeling is therefore needed in order to improve the accuracy of the predictions through the calibration of important physiological parameters. The Bayesian statistical framework to inverse problems is a powerful approach that relies on posterior sampling techniques, such as Markov chain Monte Carlo algorithms. The generation of samples re-quires many evaluations of the cardiovascular parameter-to-observable model. In practice, the use of a full cardiovascular numerical model is prohibitively expensive and a computational strategy based on approximations of the system response, or surrogate models, is needed to perform the data as-similation. As the support of the parameters distribution typically concentrates on a small fraction of the initial prior distribution, a worthy improvement consists in gradually adapting the surrogate model to minimize the approximation error for parameter values corresponding to high posterior den-sity. We introduce a novel numerical pathway to construct a series of polynomial surrogate models, by regression, using samples drawn from a sequence of distributions likely to converge to the posterior distribution. The approach yields substantial gains in efficiency and accuracy over direct prior-based surrogate models, as demonstrated via application to pulse wave velocities identification in a human lower limb arterial network.

Résumé. En s'appuyant notamment sur les progrès réalisés dans le domaine des méthodes numériques pour l'interaction fluide-structure, la modélisation du système cardiovasculaire progresse fortement. Elle est aujourd'hui en passe de proposer des outils numériques personnalisés pour l'aide au diagnostique,à la prédiction des risques,à l'intervention et au traitement clinique. Cependant la fiabilité des simulations est inévitablement impactée par les approximations liéesà une modélisation encore trop grossière ou partielle. L'intégration directe de données spécifiques au patient est donc souhaitableà l'amélioration de la précision des prédictions par le biais de l'inférence des paramètres physiologiques importants. Le cadre statistique du formalisme bayesien se prête naturellementà la resolution des problèmes inverses. Il s'appuie sur des techniques d'échantillonnages a posteriori, telles que les méthodes de Monte Carlo par chaines de Markov. La génération de la chaine requiert de multiplesévaluations du modèle cardiovasculaire reliant les paramètres aux l'observables. En pratique, le recourtà un modèle cardiovasculaire complet reste trop couteux. Dans ce cas, le choix d'une stratégie de calcul reposant sur des approximations par un modèle de substitution de la réponse du système rend l'assimilation des données plus abordable. De plus, comme le support de la distribution des paramètres Introduction Inverse problems, encountered in the parametric identification, data assimilation or shape optimization are ubiquitous in life sciences and related interdisciplinary fields such as biomedical, biomechanics and bioengineering applications. Being generally ill-posed, inverse identification often entails very large computational efforts as it requires to iteratively solve the equations characterizing the biological system (i.e. the forward problem) numerous times; the idea behind the process being to adjust some parameters in order to lower the discrepancy between the numerical prediction and some indirect clinical observations of the system. The computational burden of repeatedly solving the forward problem in order to get some outputs of interest is further amplified because the model (e.g. patient-specific cardiovascular model) is often nonlinear and brings in unknown high-dimensional parameter spaces, and potentially large but noisy datasets. Model inversion in the presence of measurement errors must typically take advantage of some -type of regularization (e.g., Tikhonov regularization) in order to recover the existence and uniqueness of solutions or -robust optimization method [13]. A potentially more natural setting is the Bayesian statistics.
Bayesian statistics can be comprehended as a systematic use of probability to decision making in the face of uncertainty. It is a completely probabilistic approach to inference. Probability models are set up for the unknowns and the data. The inference is made from the conditional distribution of the unknowns given the data: the so-called posterior distribution. Specification of a full probability model involves the determination of the likelihood function as well as the specification of a prior distribution which expresses probabilistically what is known about the unknowns before observing the data. The Bayesian framework has been successful in numerous fields of life sciences such as biomedical applications, biostatistics, biomedical imaging, neurosciences/brain, gene expression and proteomics, clinical trials, epidemiology and biomechanics modeling. Nevertheless, very few developments associated with cardiovascular mathematics and modeling are reported in the literature.
The accuracy of patient-specific cardiovascular simulation tools together with medical images and segmentation algorithms, has significantly increased in recent years, allowing for tremendous progress in the understanding and treatment of various medical conditions. One example of such achievements is the simulation of the full-scale fluid-structure interaction between the blood flow and the arterial walls of the systemic circulation network, thanks to the coupling of Navier-Stokes equations with elastic parietal models [49]. Nevertheless, reliability and sensitivity of this type of numerical predictions with respect to the model bio-physiological and biomechanical parametric uncertainties and errors are seldom reported. There exists for instance numerous uncertainties related to geometry (e.g., arterial wall thickness, lumen diameter in the case of stenosis model, lesion length), material properties (e.g., permeability, elasticity, compliance or Young's modulus), blood rheology/viscosity, micro-vasculature peripheral resistance (e.g., the boundary conditions), etc. . . associated with cardiovascular simulations. When available, patient-specific data are often scarce and indirect observations of the quantity of interest (the model parameters) and are altered by measurements noise and/or averaging acquisition procedures (e.g., clinical and medical imaging). Otherwise, population-averaged values are commonly used in place of patient-specific data. In addition, accurate full-scale direct numerical simulations remain generally out of reach and are not of effective use in support to the clinician's diagnostics and interventions. It is then a common usage to rely on reduced-order approximations (ROM), which are computationally more efficient but also more prone to model errors. After parametric calibration, the carefully crafted low-order model requires validation. In the case of vascular hemodynamics simulations, some exhaustive comparison between 1D and 3D hemodynamics formulations have for instance been published, e.g. [8,30]. Several studies have shown the ability of the 1D formulation to capture the main features of pressure, flow and area waveforms in large human arteries using in vivo measurements [9,31,48], or in vitro experiments or 3D numerical data [8,30]. Nevertheless, deterministic ROM constructed with optimized parameters value thanks to in vivo measurements or 3D numerical data often have limited value in predicting unobserved scenarios. Moreover, they do not provide an easy access to global parametric sensitivity analyses, that are essential to guide the modeling effort.
Several works, e.g., [10,12,[36][37][38], have demonstrated the interest of incorporating generic inter-and intra-patient variability in the form of uncertainty in the modeling of the cardiovascular system, since many aleatory and epistemic uncertainties remain due to its complexity, diversity and multiscale nature. These studies have shown how to propagate parametric uncertainties into the model in order to determine confidence intervals and statistics on the simulation predictions. When the level of uncertainty is overwhelming and one holds clinical data/measurements for a given patient, the uncertainty measure of his physiological parameters may be narrowed. Provided a prior distribution of the parameters to be calibrated is available, the Bayesian inference (BI) is a rigorous approach to this estimation problem and delivers, after inference, a complete probabilistic characterization of the parameters. Within this approach, the parameters are considered as random quantities which can be sampled according to their posterior distribution using, in particular, a Markov-Chain Monte-Carlo (MCMC) sampler. This technique is straightforward to implement, but it becomes computationally intractable when the forward model (FM) is expensive to solve. A simple way of alleviating this problem is to substitute the FM with a cheaper surrogate model (SM), see e.g. [29]. One of the main limitations of the surrogate acceleration for BI is that the construction of the surrogate usually aims at minimizing the approximation error for the whole prior distribution, and as such is blind to the data. In cardiovascular simulations, the nonlinearity of the FM and the dimensionality and range of the parameters make very costly the construction of an accurate SM over the whole prior distribution. When the observations are informative enough, the posterior distribution concentrates on a small fraction of the prior support and the SM needs to be accurate essentially only over the support of the posterior, as low posterior areas are not interesting. Local adaptive surrogates [27] are good candidates for the purpose of accelerating MCMC samplers, but these approximations are "memoryless" and the resulting sampling strategy may lack performance. The present work addresses some of these issues and proposes a novel formulation designed to accurately approximate the FM over the whole support of the posterior distribution of the parameters. A sequence of global, data-driven, polynomial SM is derived and used for the MCMC sampling of the posterior. Once the SM converges, in the sense of the measure of the (approximated) posterior, no subsequent call to the full FM is required. The proposed inverse modeling approach focuses on the efficient characterization of sensitive reduced-order model parameters based on clinical data. In addition to the obvious asset of uncertainty reduction, thanks to the calibration, the stochastic modeling provides patient-specific quantitative information about the confidence in simulation predictions of unobserved cardiovascular states and points to the most influential parameters, thanks to calibrated global sensitivity analysis.

Parametric Bayesian inference based on surrogate reduced-order models
Here, we introduce our notations and recall the basics of Bayesian inference for parameter identification. Then we show how the formulation may be loosened to speed up the inference. We consider a forward model that maps some unknown parameters θ consider random to some observations, or the quantities of interest, z derived from the forward model solution y. The vector of random parameters to be inferred is written as θ ≡ θ(ω) = (θ (1) , . . . , θ (n θ ) ) ∈ I θ ⊆ R n θ , The FM solution (considered discrete) is a n θ −variate functional the parameters, y : I θ → I y ⊆ R ny , so that the predicted quantities of interest are also random functionals: z : I y ⊆ R ny → I z ⊆ R n d . Denoting d ∈ R n d the dataset of observations (or measurements), the Bayes' formula writes: where prior information about θ is encoded in the prior density π prior (θ), π (d|θ) is the likelihood function that incorporates both the data and the forward model and π post (θ|d) the sought parameters posterior density. The likelihood results from the combination of the costly deterministic forward model M : R n d × R nx → R ny , an observation operator: G : R ny → R n d and statistical models for measurement noise and model error. Assuming additive measurement noise, mutually independent from the parameters, we have: where the distribution π ε of ε is prescribed. In this case, we will assume that the noise model involves no hyperparameters, so the likelihood function becomes: The likelihood function therefore contains a stochastic source term whose representation must encompass the response of the deterministic forward model over the support of π prior (θ). In practice, no closed form analytical expression exists. Posterior moments, expectations or maximum a posteriori values must be estimated via sampling methods such as MCMC, which require many evaluations of M(θ, x) and can range from a few thousands to a few millions depending on the problem. One obvious way to accelerate this computation is to substitute a faster and cheaper reduced-order model (ROM) M rom to the full deterministic forward model. Here we refer to numerical reductions with respect to the x quantity, which denotes the variables or parameters of the system whom quantified representation is known (e.g. spatial coordinates). Previous works have relied on ROMs to solve more efficiently computational inverse problems, exploiting either the deterministic, the frequentist [20], or the Bayesian approach [16]. For instance, deterministic reductions can be introduced using a Proper Orthogonal Decomposition [41] or a Reduced Basis method [32], dimension reduction, simpler geometry, etc...). A complementary approach consists in introducing a stochastic approximation that will alleviate the systematic sampling of the (costly) forward model, e.g. [29]. In our case, we build a stochastic approximation of the ROM, i.e. a surrogate model notedM rom , this time approximated in terms of its dependence to the uncertain parameters θ. In this setup, the surrogate-based posterior distribution now takes the following form:π The surrogate model should approximate the quantity of interest as accurately as possible. Classically, it is constructed from a finite set of model predictions for selected parameter values in a "off-line" stage, and subsequently used "in-line" during the sampling stage. The construction and querying of the surrogate model should be computationally efficient with sufficient predictive capabilities. Different methods are available to construct such surrogate models; examples are Gaussian processes [33], support vector machines [43], stochastic interpolation [44], or stochastic spectral methods such as polynomial-based representations. The choice of the surrogate model should be made to guarantee the accuracy of the approximate posterior distribution incurring from the substitution of the model prediction in Equation (4). Previous works have shown that the approximate posterior distribution error (in Kullback-Liebler divergence norm) is bounded by the surrogate model error measured in the mean square sense induced by the prior measure [7,29]. Therefore, polynomial chaos surrogate models [17,25,50] are well suited as they indeed minimize the error norm z −ẑ 2(πprior) = E πprior {(z −ẑ) 2 } 1 2 . In this paper we will focus also on the use of polynomial expansions, that are well-suited for smooth dependences, constructed by discrete least-squares type approach. The surrogate model then consists in a linear combination of prescribed multivariate polynomials in θ. We shall assume that z(θ) is a second order random vector for θ ∼ π prior (θ); its polynomial surrogate will write: where a γ are the unknown expansion coefficients, Λ p is an multi-index set (to be defined) for multi-index γ = (γ 1 , . . . , γ n θ ) ∈ N n θ . The polynomial approximation space is P Λp ≡ span{ψ γ | γ ∈ Λ p }, where the multivariate polynomials are often constructed as product of one uni-variate polynomials, ψ γ (θ) = having degree γ i . We restrict ourselves to isotropic tensor-product spaces total degree (TD) p leading to Λ TD p = {γ ∈ N n θ : ||γ|| 1 ≤ p}. We shall assume that these polynomials set forms a basis for any subsequent measures associated to θ, in particular the successive approximations of the posterior density. For simplicity in the construction of the original polynomial basis, we assume prior independence of the prior parameters. Posterior parameters however may be dependent or at least correlated. In this latter case, a numerical conditioning of the parameters samples (e.g. Cholesky-type decomposition), mapping the parameters into centered and normalized uncorrelated coordinates will facilitate the polynomial regression.

Adaptive weighted regression construction
The construction of accurate polynomial surrogates of general functions over large (prior) supports remains very challenging and demanding. Moreover, it is very inefficient in the context of Bayesian inference when the data are informative, i.e. when the posterior density of the inferred parameters highly concentrates from the initial prior density. A more efficient and accurate approximation can be obtained by constructing a "posterioradapted" surrogate model that would minimize the z −ẑ 2(πpost) , that is the error norm induced by the posterior. A legitimate concern is the one of the construction of a surrogate in the important region of the posterior distribution before actually characterizing the posterior. As this measure is not known from start, we propose a sequential strategy where a sequence of surrogatesM Several approaches are available for building polynomial approximations: e.g. interpolation, projection, regression, compressive sensing, etc... that only require observations of the mapping θ → z(θ). All of theses approaches have a computational cost that is dominated by the FM evaluation. Therefore, the goal is to limit the total number N of discrete samples of the mapping needed to construct the surrogate model. Because the parametric inference is driven by random sampling, we choose to rely on iterated polynomial regressions over sets of randomly generated points. Moreover, we make the choice of numerically emphasizing the samples for which the current surrogateẑ (k) predicts well z. This is ensured by associating scalar weights ρ j > 0 to the samples in the resolution of polynomial regression problem. In the following, we explain how the samples are incorporated, weighted and used in the construction of the surrogates. Let us assume, for the sake of comprehension, that we dispose of a current surrogate polynomial modelM (k) rom at some step (k) of our sequential process. This surrogate, constructed on the basis of a set of previous weighted model simulations S (k) , is fast to evaluate and serves the purpose of the MCMC sampler that produces a new set of parameters drawn fromπ (k) post (θ|d) according to the Bayes' rule of Equation (4). One may then consider a few number n (k) of well-mixed parameters samples θ (k),1≤j≤n k of this chain, and performs the corresponding simulations using the ROM FM to compute z(θ (k) , j). For each sample, one can then quantify the current surrogate error using a discrepancy measure ∆ There exists different ways of deriving sample weights from the discrepancy measure, but the main idea is to make it inversely proportional with some numerical bounding capability, i.e. ρ j ∼ max( , ∆ j ) −1 [23]. If all ρ j are large, it means that the surrogate is quite accurate in the sampled approximate posterior region, implying in turn thatπ (k) post approximates well π post . On the contrary, when ρ j is low, it indicates that the sample θ (k),j is not a reliable observation point as it is drawn from an inaccurate posterior. The sample weight completes the observation triplets constituting the set ..,n (k) of observations collected at iteration k. This set is added to the current set A of all observations collected so far: At this stage, a choice needs to be made in order to select among this growing database A, the subset of samples that will be used to construct the next polynomial surrogateM (k+1) rom at step (k + 1). Again, there is no unique bandwidth selection choice, and one may only consider selecting among the more recently generated samples, select samples with larger trust indexes across multiple previous iterations, . . . Irrespective of the selection strategy, let us denote S (k+1) ≡ {(z j , θ j , ρ j )} j∈J (k+1) ⊆J the selected subset of A. Depending on the number and quality of the samples in S (k+1) , the approximation space must be tailored. A simple option consists in relating the total order p of the polynomial approximation to the number of elements in S (k+1) .
In this work we used the following strategy. First, to construct the initial (k = 0) surrogate model we set p (0) = 1 and generate n 0 samples θ j drawn directly from the prior π prior (θ) to obtained the initial set of triplets using arbitrarily ρ j = 1. The size of the initial set is n 0 = 3 × n pol (p (0) ) where n pol (p) is the dimension of the TD polynomial space at degree p. The whole set is used to compute the initial surrogate, S (0) = A. Then, at each iteration (k), we generate n f × n pol (p (k) ) new samples from the current surrogate posterior, for some n f ≥ 1 integer. The polynomial degree of the next surrogate is selecting according to the rule For the construction, we attempt to use asymptotically the last half of the collected samples in A using the definition J (k+1) = {n − max(n/2, n f × n pol (p (k+1) ), . . . , n}. Here, we choose n f = 2 to ensure at least twice as many samples as polynomial coefficients. The value of n f can be increased to improve the stability of the regression in particular when p (k+1) becomes large. The surrogate modelM is therefore obtained by solving a weighted least squared residual minimization problem: The surrogate model is then incrementally optimized to minimize the response error in the posterior norm. If z presents sufficient regularity with respect to the parameters, one may expect fast convergence with a reasonably low polynomial degree. This translates into significant reductions in the number of FM evaluations. The convergence may be monitored numerically, for instance observing the samples weight distribution. Finally, the iterations are continued until convergence is reached or a maximum number of iterations (k max ) or simulations |A| > N is attained.

Pulse wave velocities identification in cardiovascular systemic arterial networks
In the following, we propose to apply our technique to the characterization of a subject-specific hemodynamic model of the pulse wave propagation in a vascular network. Quantifying the relationship between the physical properties of the cardiovascular system and the shape of the arterial pulse wave is important because the latter carries valuable information for the diagnosis and treatment of cardiovascular diseases. The coupled interaction of the blood flow with the network of compliant arteries is a (nonlinear) fluid-structure interaction (FSI) problem with a closed distribution network, pulsatile fluid flow with complex internal dynamics and multiple reflections leading to relatively complicated pulse wave patterns [46].

Arterial stiffness
Cardiovascular diseases, such as atherosclerosis, aneurysm, and dissection all involve significant degeneration of the arterial wall tissue, e.g., deposition of calcified materials, reduction in elastin content, plaque formation [45]. These degenerations induce changes in the arterial stiffness, i.e., the capability of the vessel to accommodate and damp the pressure waves generated by the sudden change in blood pressure due to the heart systolic contraction. Arterial stiffening plays a key role in the development of cardiovascular diseases. Even if the exact mechanism by which it leads to increase in pulse pressure and systolic hypertension is still a controversial topic [34], it remains a valid predictor of cardiovascular morbidity and mortality. The problem is those stiffness properties of the aortic walls are not well known and subject to large (spatial) variabilities among individuals [35]. In fact, pulse wave velocities not only fluctuate during a cardiac cycle but also vary along the arterial tree due to the natural stiffening of the arterial walls toward distal locations. Arterial stiffness is hard to measure in vivo [47] and require numerical substitutes [1]. Common medical practice is to infer it from an indirect measurement of pressure (averaged) pulse wave velocity in the large vessels of the arterial tree -carotid-femoral or brachial-ankle assessments are common -thanks to medical imaging (Doppler Ultrasound, MRI, CTComputed Tomography) and the Moens-Korteweg equation [22]. For instance, the ankle-brachial arterial pressures index which reflects peripheral arterial obstruction due to atherosclerosis, is an established vascular marker for the diagnosis of peripheral arterial disease which affects more than 20% of the elderly population. In this case, accurate modeling of the arterial pulse wave velocities in the limbs are critical to assessing the revascularization of the lower extremities and its impact on a major cardiovascular risk [18]. Some UQ studies exist that have confirmed that the inherent aortic stiffness uncertainty (either due to biological variability or age-related) had a strong influence in terms of standard cardiovascular outputs such as pulse pressure, brachial to radial pulse pressure amplification index and some waves transit times across the arterial tree, e.g. [10,12,24,28,37,51]. Nevertheless, in these studies, the system under consideration is sometimes too generic and the uncertainty modeling too simple and idealized, e.g. disregarding the correlations and physical constraints among parameters. Other studies have considered the assimilation of patient-specific clinical geometrical or functional data/measurements together with the use of a reduced-order model of the arterial circulation by relying on a deterministic or statistical formulation [3-6, 11, 13, 26]. When measurement error is included in the observations and uncertainty considered in the forward model describing the effect of the sought-after parameters onto the outputs, it results in a stochastic inverse problem whose solution is a probabilistic description of the parameters [3], in contrast with a point-wise estimate as in inverse problems. We propose to apply our approach to the calibration of an arterial hemodynamic model. We treat pulse wave velocities and arterial lumen area of the arteries at rest as piecewise-constant but unknown quantities modeled a priori as independent random variables (one for each artery component considered among {A i } i=1...d=7 arteries). With the help of a simplified hemodynamic model together with patient-specific non-invasive local measurements of the vessel motion and blood flow velocity, we then apply our iterative approach in a Bayesian framework in order to efficiently calibrate these pulse wave velocities.

Governing equations
We consider a portion of the systemic arterial tree containing a finite number of arteries: we have a network of thin, deformable, and axisymmetric arterial segments filled with blood, taken as an incompressible Newtonian fluid [40]. The formulation, for each arterial segment, based on the conservation of mass and momentum laws and on Young-Laplace equation, is: where t denotes time, x ∈ D is the axial coordinate along the arterial centerline, A(x, t) is the circular crosssectional area of the lumen, u(x, t) and p(x, t) are average velocity and internal pressure, respectively, ρ is blood density. The term f is the friction force per unit length and is related to the velocity profile through f = −2µπ α α − 1 u, where µ is blood dynamic viscosity and α ∈ [0, 1[ is a correction factor accounting for the nonlinear integration of radial velocities in each cross-section [42]. The underscript 0 denotes quantities at rest. The β parameter is a measure of the arterial wall stiffness related to its mechanical behavior: β = √ π h 0 E/(1 − ν 2 )A 0 , where h 0 is the reference arterial wall thickness, E is the Young's modulus and ν = 1/2 is the constant Poisson's ratio. It is also related to the pulse wave velocity c through the following relation: Note that c and c 0 increase with increasing elastic modulus and wall thickness, and decreasing luminal area. At some point in the following, due to the lack of knowledge of the pulse wave velocity at rest c 0 , it will be modeled as a random quantity which implicitly denotes the uncertainty in Young's modulus. We refer the reader to [15] for more details about the mathematical form of this system of equations.

Numerical solver
In this work, we rely on a discontinuous Galerkin (DG) method with a Legendre spectral/hp spatial discretization [21] and a Riemann solver of Roe for the calculation of the upwind fluxes, to solve the system of Equations (7). It is a very efficient scheme for high-order discretization of convection-dominated flows as it propagates waves at different frequencies without suffering from the excessive dispersion and diffusion errors.
Detailed mathematical formulation and analysis about the hyperbolic system and the numerical scheme may be found elsewhere, e.g. [9].
We use an explicit second-order Adams-Bashforth scheme with a time step dictated by the CFL condition [21]. In the case of polynomial approximation of order p = 4 within each cell, pulse pressure convergence is reached for an average mesh resolution of 20 cells per unit meter. Dirichlet-type time-dependent boundary conditions at the inlet of the domain are imposed thanks to subject-specific echo-tracking clinical data. At each exit of the simulated arterial tree, a simple terminal reflection is imposed via a scalar coefficient R ∈ [0, 1] that relates the backward to the forward characteristic information. It enables the simulation of reflected waves induced by the resistances of the missing peripheral arterioles and capillaries. This coefficient may be interpreted as an avatar of a lumped parameter model made of a single resistance [2,13].  Figure 1. Network of the seven main human lower limb arteries. Resistance symbols represent partially reflecting outflow boundary conditions. Some medical imaging data relative to local temporal evolutions of cross-sectional area and blood flow velocity are available in the arteries in blue; only flow velocity data are available for the artery in red; nothing is measured in the arteries in grey color.

Numerical results
In the following, we are interested in a subject-specific data assimilation on an arterial human left lower limb model. Here we consider the same seven-artery simple network model described in details in reference [13], where -only the largest arteries are retained and physiological properties at rest are considered piecewise constant, see Figure (1) and for some of which -we hold some medical imaging data. The non-invasive measurements are collected on the male subject thanks to ultrasound echo-tracking technique [14]. The data are cross-sectional lumen area A(x center , t) changes measured at the approximate center of the iliac, femoral and poplital arteries {A i } i=1,2,4 , and blood velocity changes u(x center , t) measured at the same locations and at the center of the anterior tibial artery {A 6 } as well, see Figures (1,5). These data are collected non-simultaneously for the same subject during two to three cardiac cycles at a time. They exhibit variabilities both due to the imaging acquisition tool and practitioner technique varying accuracy and the subject physiological changes (e.g. heart beat) during the clinical examination. In order to reduce these uncertainties and also to make the application more challenging, we condense the temporal information into scalar quantities. Contrarily to the study of El Bouti [13], where entire time signals were put to use, here the data are utilized in the form of scalar relative arterial changes measured over time: d ≡ {d l } l=1:n d ∈ R n d =7 with: where φ {A i(l) } is the l −th measured quantity at the center of the artery {A i(l) }, and t sys and t dia relate to the cardiac systolic (diastolic) time at which the measured quantity reaches its maximum (minimum) respectively. For the data we have access to, the measured quantity is either the lumen area or the blood flow velocity. Due to the lack of information, measurements statistical description is evaluated from the collected time signals and modeled as = N (µ d , Σ d ). In particular, Σ d is directly approximated from the signals temporal deviation which is a very crude assumption that does not account for instance for measurements errors due to approximate spatial acquisition sites and non-simultaneous acquisition times. The statistical description is summarized in Table 3. The chosen values are such that the coefficient of variations of these measurements exhibit a 40% dispersion. The main parameters to be inferred are the pulse wave velocities at rest c 0 ∈ R d=7 in all arteries {A i } i=1,...,d=7 of the network. Meanwhile, we also infer the resistance parameters at each outflow boundary condition, i.e. R ∈ R d=4 in arteries {A i } i=3,5,6,7 and the lumen cross-section of the arteries at rest that have not been measured, i.e. A 0 ∈ R d=4 in arteries {A i } i=3,5,6,7 . All the parameters, θ ∈ R d=15 , are expressed a priori as independent normal random variables with p(θ) ∼ N (µ θ , Σ θ ). Mean values µ θ and standard deviation Σ θ are chosen from the literature [39], and correspond to coefficient of variations of the order of c v θ ∼ 15%. In practice, these statistics insure positivity and satisfy the hyperbolicity condition of the system [51] almost surely. The iterative surrogate model of the system response is adaptively constructed based on the following numerical tools: -a deterministic parallelized DG 1D fluid-structure interaction hemodynamic solver (here with a typical resolution: h ∼ 0.2 cell/cm; and a Legendre polynomial approximation basis of order p = 4 in each mesh cell); -an automated data samples selector and preconditioner; -a set of orthonormal polynomial (here Hermite) libraries of arbitrary degree p and cardinality n pol (here relying on tensor-products constructed to satisfy a total degree expansion) and -a weighted least-square w−LS solver; the posterior sampling of the parameters is handled thanks to parallel MCMC chains (Metropolis-Hastings scheme) of 3.2 · 10 5 samples.
For this application, the forward model is too complex to make possible the existence of an exact solution mapping the quantities of interest to the uncertain inputs. Due to this lack of reference solution, the convergence of the numerical approach may be apprehended by monitoring the evolution of the surrogate model as well as the posterior characteristics of the parameters to be inferred. In Figure 2-(a), we illustrate the convergence of the adaptive surrogate construction for N = 2500 1 and p (0) = 1, reporting the evolution of the samples trust-index value (colored dots) and their average (red solid line) over the new data completing at each iteration (k) the total data set A. The results are reported as a function of the cumulated number of exact model solves; also reported is the evolution of the surrogate model order p (k) and the corresponding set of selected J (k) samples in Figures 2-(b-c). We see from the plot that the averaged trust-index globally increases, denoting the progressive improvement of the surrogate. Vertical black dotted lines materialize the points at which the surrogate order is automatically incremented. We notice patent increases when the surrogate order is augmented. Posterior statistics convergence (not displayed here) show that mean values are nicely calibrated and reach asymptotically statistically stable values after about ∼ 500 model evaluations. Concerning the std results, the data are informative enough so that the uncertainty ranges are often reduced.  Other inferred parameters distributions are not represented as they are less influent. Globally we can say that the measurements in {A 2,4,6 } arteries are informative enough in the sense that the stiffness uncertainty of the material system is lowered except in distal arteries {A 3,5,7 }. It is particularly true for the top {A 1 } artery. Located upstream of the network and contiguous of the deterministic inflow boundary condition, the system response in this artery is not prone to large uncertainties. Feeding on the downstream measurements, the arterial stiffness of this vessel is therefore sharply calibrated. This directly benefits the calibration of the connected downstream main arteries, in particular, {A 2,4,6 }. In fact, we observe very noticeable correlations between the stiffness of arteries {A 1,2 }, {A 1,4 } and also {A 2,4 }. This makes sense due to the topology of the network, see Figure 1. We also notice that joint measures between the stiffness of {A 6 } and each of {A 1,2,4 } get nicely informed from the bayesian updating. The stiffness of artery {A 7 } is somewhat updated but one    Figure 4 are an example of the capability of the method to gradually reach over remote posterior regions of interest. Once the right region is reached, more available samples render possible the construction of a higher-order polynomial approximation of the response, which makes the surrogate-based posterior sampling more efficient. The Bayesian formulation also gives access to a valuable point-wise estimate of the inferred parameters: the maximum a posteriori probability (MAP) estimate that is defined as the mode of the posterior distribution: θ MAP = argmax θπpost (θ|d). In the following, we compare the results of the numerical model relying on the MAP hemodynamic parameters with those calibrated with a different approach. Another interesting subsequent study would have been to perform an uncertainty quantification of the effect of the entire probabilistic description of the calibrated parameters onto the quantity of interest. Finally, the adaptively weighted regression approach yielded substantial gains in efficiency and accuracy over methods using direct prior-based surrogate models. A practical comparison -not presented here -was made for a similar but simpler version of the lower limb system introduced previously, i.e. with only seven stiffness parameters to calibrate. Our adaptive surrogate model approach was about one order of magnitude more efficient in terms of number of calls to the deterministic solver compared to a calibration based on a global surrogate model constructed over the prior parametric domain from a Gauss-Hermite Smolyak-designed sparse grid.

Comparison with an evolutionary optimization approach
A different pathway to the calibration of hemodynamic model parameters has been used in a previous work [13]. In that study, a Covariance Matrix Adaptation Evolution Strategy (CMA-ES) had been chosen. Is it an efficient stochastic, derivative-free evolution strategy with an auto-adaptive covariance matrix for numerical optimization of non-linear or non-convex continuous optimization problems, based on the principle of biological evolution [19]. Table 2 summarizes the different calibrated parametric results and compares them with nominal values. We mention that the CMA-ES-optimized parameters are averaged values obtained over ten reiterated optimization procedures. The first column designating the artery number also indicates what type of subjectspecific data was measured for each artery (either cross-sectional area and/or blood velocity). We see that the MAP parameters obtained from the Bayesian inference are close to the ones obtained via the optimization procedure, especially in the upper part of the network where data are available. Discrepancies do occur in the lowest section, i.e. the tibial arteries. These differences may be interpreted in the light of several factors. First, it should be noted that optimized parameters in this section are subject to larger variabilities due to the lack of measured data [13]. Moreover, one should keep in mind that the information content of the data provided to the statistical inference is more limited. Indeed, it only provides relative magnitude scales with no information about phasing nor timing. Finally, figure 5 compares the predictions based on the nominal and MAP parameters values with the subjectspecific measurements. The improvement is patent with the calibrated parameters, except for the cross-sectional lumen area in the iliac and the tibial arteries. Indeed, in the first case, the imposition of the upstream boundary condition data in the form of the measured area, is sufficient to drive the dilatation of the iliac artery to the right profile, even with parameters nominal values. In the second case, the measurements are not informative enough to induce real changes in the lumen area profile of the tibial arteries. Predictions obtained from the optimization procedure, see Fig. 9 in [13] look very similar, in particular for the first, second and third generations of arteries of the bifurcating network.

Conclusion
We have developed a computational approach in which computationally intensive cardiovascular models with unknown physiological parameters can be approximated and then calibrated via sample-based Bayesian tools and subject-specific data. To make the procedure computationally tractable, the full forward cardiovascular model or parameter-to-observable map are approximated at the deterministic and stochastic levels. More specifically, in this work, we introduce a novel adaptive pathway to construct a series of parametric polynomial surrogate models gradually adapting to minimize the approximation with respect to the posterior measure of the parameters. In our experiments, the approach yields substantial gains in efficiency and accuracy over prior-based surrogate models. The approach is applied to subject-specific pulse wave velocities identification in a human lower limb arterial network. Pulse propagation in a vascular network is a challenging multi-physics problem dominated by numerous high-speed pressure waves traveling forward and backward in the network, due to multiple bifurcations and terminal reflections. Despite this complexity, the numerical method succeeds in calibrating the physiological parameters adequately and the most probable parameter values compare well with the ones obtained using an alternative approach based on an evolutionary optimization algorithm. In addition, our approach provides richer information in the form of a complete statistical characterization of the inferred parameters. This result is remarkable considering that: a) the deterministic reduced-order model is quite simplified compared to a more sophisticated full three-dimensional fluid-structure interaction model, both in terms of network geometry, scales, and numerical modeling assumptions; b) the subject-specific real clinical data of different nature are scarce and incomplete; and c) the computational budget is relatively modest regarding the high-dimensional parametric space to explore.