Roch Molléro

Xavier Pennec

Hervé Delingette

Nicholas Ayache

Maxime Sermesant
email: maxime.sermesant@inria.fr

Population-Based Priors in Cardiac Model Personalisation for Consistent Parameter Estimation in Heterogeneous Databases

Keywords: Cardiac Electromechanical Modeling, Parameter Estimation, Parameter Selection, Personalised modeling

Personalised cardiac models are a virtual representation of the patient heart, with parameter values for which the simulation fits the available clinical measurements.

Models usually have a large number of parameters while the available data for a given patient is typically limited to a small set of measurements, thus the parameters cannot be estimated uniquely. This is a practical obstacle for clinical applications, where accurate parameter values can be important. Here we explore an original approach based on an algorithm called Iteratively Updated Priors (IUP), in which we perform successive personalisations of a full database through Maximum A Posteriori (MAP) estimation, where the prior probability at an iteration is set from the distribution of personalised parameters in the database at the previous iteration. At the convergence of the algorithm, estimated parameters of the population lie on a linear subspace of reduced (and possibly sufficient) dimension in which for each case of the database, there is a (possibly unique) parameter value for which the simulation fits the measurements. We first show how this property can help the modeler select a relevant parameter subspace for personalisation. In addition, since the resulting priors in this subspace represent the population statistics in this subspace, they can be used to perform consistent parameter estimation for cases where measurements are possibly different or missing in the database, which we illustrate with the personalisation of a heterogeneous database of 811 cases.

INTRODUCTION

Personalised cardiac models are of increasing interest for clinical applications [START_REF] Kayvanpour | Towards Personalized Cardiology: Multi-Scale Modeling of the Failing Heart[END_REF][START_REF] Baillargeon | The living heart project: a robust and integrative simulator for human heart function[END_REF][START_REF] Smith | euheart: personalized and integrated cardiac care using patient-specific cardiovascular modelling[END_REF] . To that end, parameter values of a cardiac model are estimated to get a personalised simulation which reproduces the available measurements for a clinical case. Then the personalised simulations can be used for advanced analysis of pathologies. In particular, recent works have been successful in predicting haemodynamic changes in cardiac resynchronization therapy [START_REF] Sermesant | Patient-Specific Electromechanical Models of the Heart for Prediction of the Acute Effects of Pacing in CRT: a First Validation[END_REF] , ventricular tachycardia inducibility and dynamics [START_REF] Chen | Biophysical modelling predicts ventricular tachycardia inducibility and circuit morphology: A combined clinical validation and computer modelling approach[END_REF] , as well as in detecting and localising infarcts [START_REF] Duchateau | Infarct localization from myocardial deformation: Prediction and uncertainty quantification by regression from a low-dimensional space[END_REF] using 3D personalised models.

Most cardiac models depend on many parameters, especially in 3D models where each parameter can take a different value at each node or region of the mesh, [START_REF] Marchesseau | Fast parameter calibration of a cardiac electromechanical model from medical images based on the unscented transform[END_REF][START_REF] Kayvanpour | Towards Personalized Cardiology: Multi-Scale Modeling of the Failing Heart[END_REF][START_REF] Balaban | High-resolution data assimilation of cardiac mechanics applied to a dyssynchronous ventricle[END_REF][START_REF] Finsberg | Estimating cardiac contraction through high resolution data assimilation of a personalized mechanical model[END_REF] . The number of parameters can be up to tens of thousands while on the other hand for an individual patient, the available clinical data is usually sparse and typically consists of a set of global measurements (such as ejection fraction, systolic and diastolic pressure) and possibly some imaging data with a considerable degree of noise or blurriness. Consequently, many sets of parameter values can lead to a simulation which can reproduce the available clinical data. The parameter estimation problem is an ill-posed inverse problem and all the parameters cannot be uniquely estimated from the measurements only.

A classic technique to estimate relevant parameters in this context is the use of prior probabilities distributions over the parameter values [START_REF] Neumann | Robust image-based estimation of cardiac tissue parameters and their uncertainty from noisy data[END_REF][START_REF] Konukoglu | Efficient probabilistic model personalization integrating uncertainty on data and parameters: Application to eikonal-diffusion models in cardiac electrophysiology[END_REF][START_REF] Brodie Aj Lawson | Unlocking data sets by calibrating populations of models to data density: A study in atrial electrophysiology[END_REF] . In this framework, parameters values come from a (usually Gaussian) prior probability distribution, which represents some knowledge or beliefs about the distribution of parameters. Then, given a set of measurements, a vector of Maximum A Posteriori (MAP) parameter values is estimated, which realises an optimal trade-off between its likelihood in the prior probability distribution and the error of fit of the simulation. The set of Maximum A Posteriori values is usually smaller than the set of values for which the simulation fits the measurements and possibly unique. A challenge to estimate relevant parameters is then to define accurately the prior probability distribution, which should be ideally as close as possible to the "true" distribution of parameters.

Here we explore an original approach which we call Iteratively Updated Priors (IUP), to adress the chicken-and-egg problem of defining an accurate "prior" when not other information is available beside the currently considered data. It consists in performing successive personalisations with priors, where the priors distribution is set from the distribution of personalised parameters in the previous personalisation. The rationale is twofold: first, if the algorithm converges, all the cases are personalised through a MAP where the prior distribution is the distribution of the population parameters itself. Intuitively this opens the possibility, for cases where the value of a parameter is unobservable, to regularise the fitting by using population information from other cases in the database where the value was observable. Secondly, when there are directions in the parameter space in which the simulated measurements do not vary (i.e. an unobservable parameter direction), the use of priors promotes the estimation of parameter values which are closer to the prior mean, which in turn makes the penalty (from the prior) stronger in this direction at the next personalisation. Intuitively, this process will tend to group parameters onto directions which are observable, and reduce the spread into directions which are not. In practice, we will see this algorithm leads the parameters to lie on a linear subspace of smaller dimension (See Figure 1 ).

In Section 2 we present the IUP algorithm in relation to two different mathematical tools. First, we express parameter estimation as a Maximum A Posteriori (MAP) estimation with prior probability distributions. Secondly, we use the Iteratively Reweighted Least Square (IRLS) algorithm from sparse regression to show that the IUP algorithm, when it converges, performs the minimisation of population-wide cost function with a sparse regulariser on the number of dimensions of the parameter space. We introduce two types of updates of the parameters of the (Gaussian) prior probability: "Full Matrix", where no assumption is made on the distribution and the covariance matrix, and "Diagonal Matrix" where the covariance matrix is supposed diagonal, which assumes independence of parameters in the distribution of the population. We demonstrate that in the first case the algorithm leads to the automatic selection of parameter directions in which parameter values are set to a constant value in the population, and in the second case to the automatic selection of parameters which are set constant.

In Section 3 we present results of the IUP algorithm on a database of 137 complete cases where the same 4 measurements are available. We perform the estimation of 6 parameters of a cardiac model likely to vary in this population, which also makes the problem ill-conditioned. We observe that the IUP algorithm leads to a trade-off between the number of dimensions of the final parameter set, and the mean error of fit of the 4 measurements across the whole database of personalised simulations. Then we show that if we impose a high goodness of fit of the personalised simulations, the algorithm leads to the selection of a linear subspace of minimal dimension in which for each case of the database there are parameter values where the simulation fits exactly the measurements. We discuss how this algorithm can be used to find a relevant subspace for personalisation from the modeling point of view. Finally we extract a relevant subspace of dimension 4, based on 5 parameters out of the original 6, in which there is a unique set of parameter values for each case of the database.

In the last section we apply the proposed framework to the consistent personalisation of a larger database of 811 cases, in which the set of measurements reported by the clinicians are heterogeneous and sometimes incomplete. Applied to this database, the algorithm leads to personalised simulations with high goodness of fit for all the cases of the population while grouping the parameters onto a parameter space of minimal dimension. Simultaneously, the resulting population-based priors computed in this parameter space lead to consistent parameter estimation for the cases where the parameters are not all observable, and are in particular influenced by their values in cases of the population where they are observable. Finally we show how external parameters such as the weight and the height of the patient can also be added to the prior distribution in order to constrain the parameter estimation with more consistency in case of unobservability.

ITERATIVELY UPDATED PRIORS AND POPULATION-BASED PRIORS

In this section we introduce the personalisation framework, based on the estimation of a Maximum A Posteriori (MAP) with prior probabilities. In order to get relevant values for each case through MAP estimation, the prior probability has to be relevant. Ideally, it should be the "true" underlying probability distribution of parameters in the population. We introduce FIGURE 1 Schematic representation (with synthetic data) of parameter estimation when the problem is ill-posed: in this example both the contractility and the stiffness are estimated from values of the stroke volume (SV). Both have an influence on the stroke volume (SV) so there are isolines of stroke volume (in grey) on which the value of the stroke volume is the same. Image a: estimation without priors, the estimated values (green) for each case can be anywhere on an isoline (grey). Image b: the estimation is performed with a (gaussian) prior (the gaussian covariance is in blue), the estimated values are grouped closer to the prior mean. Image c: the Iteratively Updated Priors (IUP) algorithm performs successive estimations where the prior is set from the distribution of estimated parameters at the previous iteration. This leads the parameters to lie on a reduced linear subspace (orange).

here the Iteratively Updated Priors (IUP) algorithm, which performs successive personalisations of the whole population, where the prior at each step is set from the distribution of the personalised parameters at the previous step.

First we present the MAP estimation of personalised parameters or a single case (Section 2.1). Then we formulate the Iteratively Updated Priors algorithm for a population of cases (Section 2.2). We then give the equations of the prior parameter update in two different cases with different assumptions on the shape of the distribution (Section 2.3) and the equations which are minimised within our optimisation framework (Section 2.4). Finally we present a link between the IUP algorithm and the Iteratively Reweighted Least Square algorithm for sparse regularisation, and show that the IUP algorithm performs the minimisation of a population-wide cost function with a penalty on the rank of the set of personalised parameters (Section 2.5). Finally we present a modification of the IUP algorithm called IUP* to impose a high goodness of fit to the model outputs to the target values for each case and at each iteration (Section 2.6).

Maximum A Posteriori Estimation of Personalised Parameters

We consider a cardiac model and a set of simulated quantities called the outputs O, such as the ejected volume and the mean ventricular pressure. We then consider a subset P M of varying parameters (such as the contractility, the aortic resistance) of the model, while the other model parameters P ′ M are supposed fixed. Given a vector of values of the parameters P M , we note O M ( ) the values of the outputs O of the model M.

With these notations, personalisation consists in estimating parameter values for which the outputs values O M ( ) are consistent with some observed values Ô of these outputs. This is an inverse problem, which can be tackled by different methods (see the review of Chabiniok et al. [START_REF] Chabiniok | Multiphysics and multiscale modelling, data-model fusion and integration of organ physiology in the clinic: ventricular cardiac mechanics[END_REF] ).

A classic approach is to formulate the problem as the estimation of a Maximum A Posteriori with prior probabilities on the parameter values. With Bayes' theorem, the posterior probability ( | Ô) of parameter values considering the observed output values Ô is proportional to the product of the conditional likelihood ( Ô| ) and the prior probability ( ):

( | Ô) ∝ ( Ô| ) ( ) (1)
In this work, we use Gaussian distributions for both the conditional likelihood and the prior probability:

( Ô| ) ∝ -1 2 ( ( )-̂ ) Σ -1 ( ( )-̂ ) ( ) ∝ -1 2 ( -) Δ -1 ( -) (2)
which correspond to the assumptions that first, the observed values ̂ are the sum of the outputs values ( ) and a Gaussian noise with covariance Σ. Second, the underlying probablity distribution of the model parameters is a Gaussian with mean and covariance Δ. Finding a Maximum A Posteriori (MAP) consists in finding (one of) the parameter values which has the highest posterior probability:

= ( | ̂ ) (3)
This can also be interpreted as finding a maximum of the joint probability ( , ̂ ) of x and ̂ with a uniform prior P( ̂ ) on ̂ , thanks to the formula:

( , ̂ ) = ( | ̂ ) ( ̂ ) (4)

Iteratively Updated Priors (IUP)

In the following we consider a population of = 1.. cases to be personalised, we note Ô = ( ̂ 1 , ..., ̂ ) and X = ( 1 , ..., ). We introduce the algorithm called Iteratively Updated Priors (IUP) which consists in successively estimating personalised parameters with the MAP estimation, and re-estimating the prior parameters ( , Δ) as the maximum likelihood of these parameters considering the estimated personalised parameters :

• Step 1: Personalisation For each case i, find = ( , ̂ | , Δ
). This consists in performing the MAP estimation in Equation 3 for each case with current and Δ, which is the estimation of personalised parameters with a given prior.

• Step 2: Re-estimation of the Prior Find ( , Δ) = ( ,Δ) (X | , Δ). This consists in finding the most likely values (maximum likelihood estimation) of and Δ given a set of personalised parameters X .

When the algorithm converges, the prior probability distribution with the resulting parameters ( * , Δ * ) is called a population-based prior based on the ̂ .

Explicit update formulas for and Δ in two cases.

We explicit here the update in the Step 2 of the IUP algorithm in which and Δ are set to the solutions of ( , Δ) = ( ,Δ) (X | , Δ). We consider two formulations of the Gaussian prior probability:

• First, a case (Full Matrix) where no assumptions are made on the covariance matrix Δ of the Gaussian distribution.

• Second, a case (Diagonal Matrix) where we suppose that Δ is diagonal. This is equivalent to the (simplifying) assumption that parameter values are independent from each other in the population.

The solutions of the maximum likelihood estimation of and Δ in the two cases are1 :

• Update of : In both cases, is set to the mean of personalised parameters :

= 1 ∑ =1 -1 (5)
• "Full Matrix" update of Δ: In this case the minimisation is reached when Δ is the covariance matrix of the personalised parameters:

Δ = 1 ∑ =1 ( -1 -)( -1 -) (6)
• "Diagonal Matrix" update of Δ: In this case the minimisation is reached when Δ is the diagonal part of the covariance matrix of the personalised parameters:

Δ = ( 1 ∑ =1 ( -1 -)( -1 -) ) (7)

Practical Implementation of the MAP Estimation.

In the following, we use a constant Σ = (N ) for the noise model in 2, where N is a normalisation vector whose coefficients are explicited later.

We find the MAP of Equation 3 through the minimization of the negative log-likelihood -2 ( ( | ̂ ), , Δ), a regularised cost-function which we call ̂ :

̂ ( , ̂ , , Δ) = ( , ̂ ) + ( ) (8)
with

( ) = ( -) Δ -1 ( -) ( , ̂ ) = ||( ( ) -̂ ) ⊘ N || 2 ( 9
)
where ⊘ is the Hadamard (coordinate-by-coordinate) division. We call ( , ̂ ) the data-fit term and ( ) the regularisation term.

We perform optimisation of this function with a derivativefree algorithm called CMA-ES [START_REF] Hansen | The cma evolution strategy: a comparing review[END_REF] , which stands for Covariance Matrix Adaptation Evolution Strategy. It asks at each iteration for the scores of points (a generation), drawn from a multivariate distribution with covariance and mean . Then, it combines Bayesian principles of maximum likelihood with natural gradient descent on the ranks of the points scores in the generation to update both and .

For each MAP estimation, we perform from 50 to (maximum) 250 iterations of CMA-ES with a population size of 20, and the algorithm also stops if all the values of √ ̂ , (which we call the data-fit term) within the search space are in an interval smaller than 0.01. The algorithm is stochastic (because the parameter values of the generation samples are drawn randomly). It iteratively enlarges its search space until it contains a minimal solution, then reduces the search space according to the best members of each generation (drawn randomly). As a consequence it does not converge to the same solution when the possible solutions are not unique. Because of this it can be used to give an (empirical) evaluation of the "non-uniqueness variability" by repeating the personalisation process multiple times.

Finally for practical reasons, since all the parameters of our model are positive, we always consider in the following the logarithm of parameter values instead of their values. This enables in particular to not have the optimisation algorithm test negative (non-physical) values in this step.

Link with Sparse Regularization

The IUP algorithm as presented in Section 2.2 has strong links with the Iteratively Reweighted Least Squares (IRLS) algorithm in sparse regression. When the IUP converges, we can show (see a more complete derivation in APPENDIX A) that it performs the minimization of a population-wide cost function, with a penalty on the rank of the set of personalised parameters:

S (X , Ô) = 1 ∑ =1 ( , ̂ ) + (10)
Where is the rank (number of non-zero eigenvalues) of Δ * , which is the covariance of the set of personalised parameters at convergence. The minimisation of thus correspond to the minimisation of the number of dimensions of the set of personalised parameters. With the Diagonal Matrix assumption, the covariance matrix eigenvectors are aligned with the coordinate axes so in this case, the IUP algorithm minimises the number of parameters which do not have a constant value in the population of personalised parameters. In order to illustrate the difference between the two types of updates, we display in Figure 2 a schematic representation of the two different behaviours.

IUP* algorithm: Equations in the Limit Case where the Noise is Considered Negligeable

The MAP estimation of Section 2.1 and the cost function minimisation in Section 2.4 makes the assumption that the measured values of an output is the sum of its true value and a Gaussian noise (Equation 2). Though this is formulation is practical and can be especially appropriate in cardiac personalisation since clinical data is often noisy, it has drawbacks. In particular, it makes the questionable assumption of having the same Gaussian noise model for all the cases in the database, whose parameter values has to be estimated or manually imposed. In addition, we might want in some context to consider the noise to be negligeable or null, and have the model outputs fit exactly the target values in the clinical data (This will be useful in Section 4).

To that end, we introduce here a variation of the IUP algorithm called IUP*, where the data-fit term is progressively reweighted during parameter estimation, until it reaches a value lower than a specific threshold. Namely, we consider a modified cost function ̂ * from the Equation 8 of Section 2.4, with an additional term :

̂ * ( , ̂ , , Δ, ) = ( , ̂ ) + ( ) (11)
In the IUP algorithm, we first initialise to 1 at the first iteration. Then for each parameter estimation of each case in the database we perform the following algorithm, which ensures that the data-fit term of every case is below 0.01:

1. optimise ̂ * ( , ̂ , , Δ, )
2. if ( , ̂ ) > 0.01 set = 2 and go to step 1.

else stop.

Then, in each subsequent iterations of the IUP algorithm, we use the final value of from the previous iteration as an initialisation. It is important to note that this formulation makes the assumption that there exists at least one set of parameter values for which ( , ̂ ) < 0.01 for which measurements are fitted very closely.

FIGURE 2

Behaviour of the algorithm for the two different updates (See legend in Figure 1 ). Image a: with the Full Matrix updates, there is no assumption on the covariance matrix so the prior distribution at the next iteration is the closest gaussian distribution (in terms of maximum likelihood) to the estimated parameters distribution. Image b: with the Diagonal Matrix updates, the prior distribution at each iteration is the closest gaussian distribution with the axes parallel to the coordinate axes.

In the final population of personalised parameters, this leads some parameters to have a constant value (here the stiffness).

From a mathematical point of view, this algorithm is equivalent to solving the following problem:

min ( ) . . ( , ̂ ) ≤ 0.01 (12)
which is a relaxed formulation of the constrained optimisation problem with the constraint that measurements are perfectly fitted:

min ( ) . . ( , ̂ ) = 0 (13)
We can also use the sparse formulation in Section 2.5 to exhibit the optimisation of a cost function over the whole database: by considering the mean of all modified cost functions over all cases, we observe that this is equivalent to solving the following problem over the complete database: (14) which is a relaxed version of the constrained optimisation problem over the whole database:

min X 1 ∑ =1 ( ) . . ∀ = 1.. , ( , ̂ ) ≤ 0.01
min X 1 ∑ =1 ( ) . . ∀ = 1.. , ( , ̂ ) = 0 (15)
At convergence of the IUP algorithm, we have ( ) = ( - * ) (Δ * ) -1 ( - * ) where * and Δ * are respectively the mean and covariance of the set of personalised parameters X (see APPENDIX A). We then have:

1 ∑ =1 ( ) = (16)
where is the rank of Δ * . Overall we are then solving:

min X . . ∀ = 1.. , ( , ̂ ) ≤ 0.01 (17)
which is a relaxed version of the constrained optimisation problem:

min X . . ∀ = 1.. , ( , ̂ ) = 0 (18)
Since at convergence Δ * is the covariance of the set of personalised parameters, in other terms, we are minimising the number of dimensions of the set of personalised parameters under the constraint that all cases have a high goodness of fit to their corresponding measurements.

RESULTS ON 137 COMPLETE CASES AND APPLICATION TO PARAMETER SELECTION.

In this section we present results on the personalisation of 137 complete cases with a cardiac model. The cardiac model we use is a fast "0D" model which we introduced in [START_REF] Molléro | Multifidelity-cma: a multifidelity approach for efficient personalisation of 3d cardiac electromechanical models[END_REF] . It is a reduced version of our in-house 3D electromechanical model, based on the implementation of the Bestel-Clement-Sorine (BCS) model 5 by Marchesseau et al. 17 , 16 in SOFA 2 . As described in [START_REF] Molléro | Multifidelity-cma: a multifidelity approach for efficient personalisation of 3d cardiac electromechanical models[END_REF] , both the 3D and 0D model share the same mechanical and haemodynamic equations, but simplifying assumptions (see APPENDIX B) are made on the geometry of the 0D model. This leads to a very fast model made of 18 equations, which can simulate 15 beats per second at a heart rate of 75 bpm. In our context it is particularly practical for experiments which involve repeated personalisations of many cases, since the personalisation of a single case with CMA-ES as described in Section 2.4 takes around 3 minutes.

We consider a set P 0D of 6 parameters of the model, and a set of 4 measurements O listed in Table 1 for which the 4 values are available for all the 137 cases.

The normalisation coefficients for this problem (in the vector N defined in Section 2.4) are 10 ml for the Maximal Volume (MV) and the Stroke Volume (SV), 200 Pa for the Diastolic Aortic Pressure (DP) and the Mean Aortic Pressure (DP). With this normalisation, an error of fit √ ( , ̂ ) lower than 1 (resp 0.1), means a personalised simulation matches volume measurements within 10 ml (resp 1 ml) and pressure measurements within 200 Pa (resp 20 Pa), which can qualitatively be considered as acceptable (very good).

Trade-off between Data-fit, Regularisation and Dimension of the Set of Personalised Parameters

We provide results of the IUP algorithm for at least 20 iterations, with both the Diagonal Matrix and the Full Matrix assumptions for values of of 0.5, 0.1, 0.05 (see names of the run in Table 2 ).

In Figure 3 , we report for each run the mean value ̄ (X , Ô) of the data-fit term ( , ̂ ) and the mean value ̄ (X ) of the regularization term ( ) at convergence of 2 www.sofa-framework.org the MAP estimations (bottom). We also perform the Singular Value Decomposition (SVD) of the set of personalised parameters and report the eigenvalues (Figure 3 , top). Then we report (rows of Table 3 ) the eigenvectors at the last iteration ( , = 1..6), and their corresponding eigenvalue ( , = 1..6).

The eigenvectors coordinates are reported as 0 in bold when the value was below 0.01. First we observe the classic phenomenon that the higher is, the higher the data-fit term (which characterizes the distance between simulated and measured values) is (Figure 3 , bottom, blue). A higher means that the noise (i.e. uncertainty on the measurements) is considered bigger in Equation 2 so estimated parameters values then tend to be closer to the prior mean than from a value which perfectly fits the measurements.

Second, that in every run of the algorithm, there is one or more eigenvalue which converges to 0 and is under 10 -4 at convergence. When the value is zero, this means that the parameters lie on a subspace of lower dimension. Here for each run, we report in Table 4 the number of eigenvectors with a eigenvalue above 10 -4 , the mean value ̄ (X , Ô) of the 4 For each run at the final iteration : number of eigenvectors for each run with a eigenvalue above 10 -4 (column ' ≥ 10 -4 '), mean data-fit value ̄ (X , Ô) across all cases and mean value ̄ (X ) of the regularization term.

data-fit term across all cases and the mean value ̄ (X ) of the regularization term.

From this table, we first observe that the mean of the regularization term is very close to where is the number of eigenvectors which do not have a value close to 0 (we use ≥ 10 -4 as the threshold for close to 0), which is consistent with Equation 10 of a cost function with a sparse regulariser on the dimensionality of the set of personalised parameters. We also observe that for the same , the data-fit term is lower for the runs with the Full Matrix updates than the Diagonal Matrix while having a lower number of non-zero eigenvalues.

More interestingly, we can observe the shape of the eigenvectors at convergence of the algorithm. In every run with Diagonal Matrix updates, we observe that the smaller the eigenvalue associated to eigenvector is, the more it is aligned to a coordinate (a vector with shape (0,...,1,...0)). This means that for these runs, the parameters associated with these coordinates have a constant values in the final set of personalised parameters. In particular, the contractility 1 has always a constant value, and either (in DIAG-IUP-05) or (in DIAG-IUP-01) is constant as well. On the other hand, there is no such phenomenon in the runs with Full Matrix.

The behaviour of the algorithm can be understood from the sparse formulation explicited in 2.5 and the cost functions which is minimised. The algorithm indeed tries to find a prior for which there is an optimal trade-off between the number of dimensions of the set of personalised parameters, and the mean value ̄ (X , Ô) of the data-fit term. Given a specific , the "cost" of lowering the dimension is an increase of in ̄ (X , Ô). With the Diagonal Matrix updates, this is done with a constraint on the prior covariance matrix to be diagonal so having a lower dimension means fixing the value of a parameter. On the other hand, with the Full Matrix updates, the algorithm can find any direction.

In particular we can observe than for = 0.1, both algorithms find a parameter subspace of dimension 4, but the data-fit ̄ (X , Ô) is higher with Diagonal Matrix updates (0.23 ≥ 2.15e-2), where both 1 and are fixed. By comparing to the DIAG-IUP-005 where only the 1 is fixed, we can interpret that the "cost" of fixing is a loss of around 0.23 on ̄ (X , Ô), which also means that the quality of some personalisations in the database is impacted.

Indeed, for the final personalisation of DIAG-IUP-01 there are at least 3 cases for which the personalisation is highly impacted ( √ ≥ 3.5) because of the fixed value of . These cases have an aortic diastolic pressure which is particularly low compared to the rest of the database ≤ 5400 . To fit this measurement, (which is the asymptotic and minimal value of the aortic pressure in the blood flow model) needs to be at least below this value (in particular in DIAG-IUP-005 and FULL-IUP-01, all measurements of these cases are almost perfectly fitted and the estimated for theses cases is ≤ 5380 ), but at this step the prior value of in DIAG-IUP-01 (which is thus the fixed constant value) is 5873 , which makes the fitting of the Diastolic Pressure impossible in these cases.

Parameter Existence and Uniqueness in the Resulting Hyperplanes

A classical question in modeling and inverse problems, is to determine which parameters are observable with respect to a specific set of measurements. For example, here we estimate 6 parameters from 4 observed outputs. If there is a linear relationship between the parameters and the outputs (i.e there exist a matrix M such as O( ) =

), the size of the kernel of the matrix is at least 2. Considering some measurements O, some parameters such as such as O = , then for any vector y in the kernel O = ( + ) as well. This means that there are at least two orthogonal parameter directions in which the parameter values are unobservable. With a non linear-model (such as cardiac models), a similar phenomenon of non-uniqueness exists locally around some personalised parameters, to the extent that the model can be approximated by its gradient. In general, there is an entire manifold of parameters for which the outputs are the same.

Using Gaussian priors on such underconstrained linear models leads a unique solution to the inverse problem, because in this case, the cost function 8 is strictly convex. To a certain extent, this can be locally true for a non-linear model and a unique specific value is promoted within the manifold of parameters for which the outputs are the same.

Finding a (possibly unique) value which is the most consistent considering a prior knowledge on the distribution of parameters is then the most interesting consideration of the MAP estimation, but this supposes to have an accurate prior. Since we are deriving the prior probability distribution from successive personalisations over the dataset itself, it is thus not possible to rely on its interpretation as a "prior knowledge" to set parameter values in the unobservable directions of this dataset.

On the other hand when no information on the statistics of the parameters is available, the only possibility is to perform a reduction of the parameter space, by forcing some parameter or parameter directions to have a fixed value. This can be done through PCA, PLS [START_REF] Rohé | Highly Reduced Model of the Cardiac Function for Fast Simulation[END_REF] or sensitivity analysis [START_REF] Marchesseau | Fast parameter calibration of a cardiac electromechanical model from medical images based on the unscented transform[END_REF] on the parameters with respect to the clinical data, and estimate coordinates of the first modes only. This also leads to a few questions regarding the reduced space of parameters, particulary on the existence and uniqueness of parameter values for which the simulation fits the target data.

Here our parameters converge into a reduced subspace in all runs of the algorithm. In particular in the FULL-IUP-005, FULL-IUP-01, DIAG-IUP-005 run, there is also a very low data-fit term across the cases. This suggests that among the undetermined directions of a parameter space, with a small enough prior, the algorithm selects an hyperplane of minimal dimension in which for each case, there is at least a set of parameter values which fits exactly the measurements.

We demonstrate this claim by analysing the resulting parameter subspaces, in these terms of parameter existence and uniqueness. To that end we perform multiple personalisations without priors (i.e. an uniform prior, which is equivalent to solving Equation 8 without regulariser, or setting = 0) where parameter are taken from within these subspaces. Because the CMA-ES algorithm is stochastic, it usually converges toward different values of the parameters at each run if there are multiple set of parameter values which minimise the cost function. We compare the following parameter subspaces:

1. The complete space H 0 of 6 parameters.

The subspace H

1 of 5 parameters of all parameters except the stiffness 1 , which is set to its final (constant) in the DIAG-IUP-005 run.

The 5 subspaces

H ( 1 , 0 ) , H ( 1 , 0 ) , H ( 1 , ) , H ( 1 , ) and
H ( 1 , ve ) of 4 parameters where both the stiffness and another parameter are set to the final prior mean value in the DIAG-IUP-005 run.

4. The 4 subspaces H 5 , H 4 , H 3 , H 2 of dimensions respectively 5,4,3 and 2, which are the hyperplanes defined by their center at the prior mean of the FULL-IUP-01, and the largest eigenvectors of the prior covariance for respectively l=5,4,3,2.

We then report both the mean error of fit √ (X , Ô) across all cases in the database, and the variability of estimated parameters across different personalisations, estimated by averaging across all cases the standard deviation of the 5 estimated values in the 5 personalisations. We first observe that the mean error of fit is 0 (lower than 10 -3 ) for only four parameter spaces: the original space H 0 , the parameter space H 1 with all the parameters except 1 , and the two parameter spaces H 4 and H 5 with respectively the 4 and 5 largest eigenvectors of the Full-IUP-01 run. This shows that they are the only subspaces which contain parameter values which fit the measurements for all cases.

Then among these subspaces, we observe that the only subspace for which both the mean error of fit and the variability of all parameters is 0 (lower than 10 -3 ) is H 4 , the hyperplane with the 4 largest eigenvectors of Full-IUP-01. In the other subspaces there is a variability from one personalisation to the other for at least the haemodynamic parameters ve , and showing that the parameters are not unique. This shows that around the prior which was found in the DIAG-IUP-005 run, 1 was the only parameter which is possible to set to a constant value without the personalisation of some cases being impacted (such as the 3 cases described in the previous section). However, once 1 is set, the variability in H 1 show that there is still an unobservable direction (especially in the haemodynamic parameters), but it is then necessarily a combination of parameters which is unobservable.

Finally, because of the possibility to find subspaces which are not necessarily aligned with the coordinates, the Full-IUP-01 run was able to find a subspace of lower dimension, 4, in which there are parameters fitting all the cases and no variability, thus reducing the space in both directions of uncertainty.

Selection of a Parameter Subspace for Personalisation

The resulting parameter subspaces in the previous section exhibit interesting properties in terms of existence and uniqueness of parameters for personalisation, so a question is how relevant are they in the context of personalisation? First, it is important to observe that there is no guarantee of uniqueness of the parameter directions which are selected by the algorithm (and the personalised parameter either). Indeed, in addition of the model being non-linear (which makes the data-fit term not convex, thus the whole Equation 8), sparse regularisations are usually not convex and have many possible solutions. In particular here, there are possibly multiple parameter subspaces on which the algorithm could converge. Visually, this can be observed in the schema of Figure 2 , where two different subspaces can be selected, each containing a unique set of parameter values for each of the 5 cases.

Secondly, there are no guarantees that the parameter spaces selected by the algorithm are relevant from a modeling or physiological point of view. In particular, when two parameters are not completely observable from some measurements in a database, their 'actual' value (if they correspond to physical parameters) in a population is likely variable. In addition, a drawback of fixing a parameter to a specific value, is that it might force other estimated parameters to vary more within the population to account for variations which would have come from the fixed parameter.

Parameter selection in the sense of setting a specific parameter direction to a constant value is then an imperfect approach in modeling, but it is unavoidable to not have variability in the estimation when no other information or statistics is available on this specific parameter direction. In this context, though the selection of relevant parameters for personalisation cannot be performed entirely automatically (because the physical meaning of parameters is ignored by the algorithms), we believe that our algorithm can help the modeler by revealing unobservable parameter and parameter directions. To that end, we recommend the following general approach, considering a set of clinical data for which observable parameters are unclear:

1. First, perform one (or multiple) runs of IUP with Diagonal Matrix updates. This is because the resulting subspace is easier to analyze, since the parameters usually have a physical meaning in the context of modeling. Then analyze the parameters which end up with a eigenvalue close to 0, and set them to a constant value if it is not incompatible with physiological considerations.

Second, perform one (or multiple) IUP with Full Matrix

updates to further select a lower dimensional subspace.

H * 0 0 0 0 0 0
We apply this approach to the current problem. First it seems physically likely that we do not have so much information on the stiffness 1 from only the 4 measurements in Table 1 , so we can reasonably decide to set its value to a constant, which we set at its value in DIAG-IUP-005. Then we perform a new run of the IUP algorithm, which we call Final-IUP-005, with Full Matrix and = 0.05 on the resulting parameters (thus in the hyperplane H 5 ). This leads to the final eigenvectors and eigenvalues:

We then select the hyperplane H * made of the 4 largest eigenvectors ( 1 , 2 , 3 , 4 ) at the end of the run, and test the existence and uniqueness of personalised parameters in this hyperplane, with the same method than in the previous section. Results are reported in Table 7 .

As expected, both the variability of parameters and the mean error of fit are 0 (lower than 10 -3 ). This means that for each case, the hyperplane contains a unique (to the extent that we can evaluate it through this algorithm) set of parameter values for which the simulation fits the measurements.

To conclude, we found a minimal parameter subspace of dimension 4, based on 5 parameters, consistent from a physiological point of view, in which to perform consistent parameter estimation. More importantly, the final parameters * and Δ * of the prior at the end of the Final-IUP-005 run respectively correspond to the mean and covariance of the estimated parameters, which we then call the population-based priors.

Finally, we point out that the resulting population-based priors will not necessarily correspond to the "real" underlying parameter distribution. Indeed, considering a population with a known parameter distribution, the priors obtained through the IUP method will correspond to this distribution only if the data enables to completely determine the distribution parameters (in which case the influence of the prior will be minimal). However, if some parameters or parameter space directions are unobservable, the algorithm will choose arbitrary directions within the space, and build priors on these directions (which will result in dimensionality reduction to a minimal subspace). It is then up to the modeller to check that the selected subspace is relevant, possibly add additional output values to personalised or, as we will see in the next section, enhance the dataset with other elements which can help determines the remaining unobservable directions.

CONSISTENT PARAMETER ESTIMATION IN A DATABASE WITH MISSING OR HETEROGENEOUS MEASUREMENTS

In this section we present the application of the proposed framework to the consistent personalisation of a large database of cases with heterogeneous (or missing) measurements. We present the personalisation of all the cases in such database with the IUP* algorithm which to thefollowing properties of the personalised simulations:

1. All the measurements of all cases are well fitted in the corresponding personalised simulations.

2. Parameters lie on a reduced subspace of minimal and sufficient dimension.

3. Parameters for cases where measurements are missing are constrained by the population-based priors in this subspace which means in particular that unobservable parameters for these cases are guided by their values in the other cases of the database where they are available.

In addition we also show in this section how the proposed framework can be extended to integrate external parameters such as the height and the weight of the patient in order to guide the estimation of parameters for cases where measurements are missing, leading to an improve consistency of parameter estimation for these cases.

Heterogeneous Database of Clinical Cases

We consider a larger database of 811 cases from different studies, hospitals and protocols. Depending on the protocol, the same measurements are not available for all the patients and even within a single study with the same protocol, some measurements can be missing. In order to accommodate the heterogeneous nature of the database, we use a heterogeneous data-fit term. Instead of using the fixed formulation ( , ̂ ) = ||( ( ) -̂ ) ⊘ N || 2 with being the vector of outputs ( , , , ) and ̂ the corresponding measurements, we build a different vector of observations for each patient. Depending on the available measurements, we use in the order of priority: As we explained in Section 3.1, the IUP algorithm with a constant noise as used in Section 3 realises a trade-off, driven by the constant , between the number of dimension of the set of personalised parameters and the mean goodness of fit over the whole database. In particular using a small , a high goodness of fit was readched for all cases, leading to the selection of a minimal and sufficient subspace in where for each case of the database there is a set of parameter values for which the simulation fits exactly the measurements.

However on this heterogeneous database, selecting a small does not lead to a high goodness of fit for all cases. As we show in the next section (Section 4.2), an IUP run on the database of 811 patients with the heterogeneous data-fit term and = 0.05 (and even as small as 0.02) leads to a parameter space of 3 directions instead of 4, and many complete cases (with 4 measurements available) are not well fitted. This is because there are cases in the database with less than 4 measurements which can be well fitted with parameters in a parameter space of dimension 3 only. Consequently the contribution to the mean error of fit of badly fitted cases is small, which leads the algorithm to remove a parameter direction which is necessary to fit other cases. A possibility would be to use an even lower , but this leads to practical numerical problems during the MAP estimation, because the value of the regularization term becomes too small in Equation 8 compared to the data-fit term. Instead, we use in the next section the IUP* version with reweighting of the data-fit term in the of the algorithm presented in Section 2.6.

Comparison between the IUP and IUP* Algorithms on the Database

Here we demonstrate differences in the behaviour of the IUP algorithm and the extended IUP* algorithm on this database. We compare two runs of these algorithms, respectively called FULL-IUP-05 and FULL-IUP * -05, based on the 5 parameters of H 1 (the stiffness is set to a constant value), a value of = 0.05 and the Full Matrix updates. We first report as in Section 3.1 the number of eigenvectors with a eigenvalue above 10 -4 at convergence, the mean data-fit value ̄ (X , Ô) across all cases and mean value ̄ (X ) of the regularization term, for each run in Table 9 .

We observe that the IUP run leads to 3 eigenvectors with an eigenvalue above 10 -4 while the IUP* run leads to 4. The mean data-fit term value ̄ (X , Ô) for the IUP* run is very low (1.35e-3) and the IUP run is also low (1.01e-2). However we

Run ≥ 10 -4 ̄ (X , Ô) ̄ (X )
FULL-IUP-05 3 1.01e-2 0.15 FULL-IUP * -05 4 1.35e-3 0.19

TABLE 9

For each run at convergence: Number of eigenvectors for each run with a eigenvalue above 10 -4 (column ' ≥ 10 -4 '), mean data-fit value ̄ (X , Ô) across all cases and mean value ̄ (X ) of the regularization term.

now show a boxplot which compares the repartition of the values of the data-fit term error of fit √ ( , ̂ ) in the population in Figure 4 : As we can see with the IUP run, there are many cases (almost one fifth of the population actually) which are not well personalised ( √ ( , ̂ ) ≥ 0.1). On the other hand in the IUP* run, by design of the algorithms all the cases have a data-fit term below 0.01.

Finally in Table 10 we also report, for the various possible values of the weight of the data-fit term (in Equation 11) at convergence, the percentage of cases in the population with these values. As we can see, this term doubled up to 4 times in some cases (thus leading to = 16).

Overall, this section shows that the extended version IUP* enables a personalisation of the complete database where all the individual cases have their measurements fitted within a specific value of the data-fit term, while parameters are gathered on a parameter subspace of minimal and sufficient dimension.

Parameter Estimation in cases of Unobservability

Here we demonstrate the impact of the population-based priors computed through the IUP* algorithm in this database on cases with unobservable parameters. To that end, we compare the estimated parameter values for cases where measurements are missing in two different parameter estimations for this database: first, the estimation of the 5 parameters without priors and second, the estimation at convergence of the FULL-IUP * -05 run of the IUP* algorithm over the complete database and the 5 parameters.

The most interesting set of values is the following: in Figure 5 , we display the (log-)values of estimated parameters as a function of Ejection Fraction (resp Stroke Volume) for the cases where only the Ejection Fraction (resp Stroke Volume) is available. We observe a classic phenomenon with priors: personalised values not only have less variance with the use of priors (blue points), but because one measurement has to be fitted, they also lie onto a space of (local) dimension 1. Indeed, this space is defined (in the case of the ejection fraction) by the equation

( ) = ( ) . . ( , Ô) = 0 (19)
It is a space of local dimension 2 (resp 3) when 2 (resp 3) measurements have to be fitted (to the extent that the model is locally approximable to its gradient). On an interpretative level, when some measurements are missing in our database and many parameter values are possible, the personalisation with priors in H * leads to the selection of a set of parameter values which maximises its likelihood in the probability distribution of the priors, or equivalently minimizes the distance ( - * ) (Δ * ) -1 ( - * ). To that extent, it performs a form of imputation in the parameter space by choosing the most likely set of parameters according to the distribution defined by * and Δ * . Since * and Δ * are the statistics of the whole population of personalised parameter values, we argue that values of unobservable parameters and parameter directions in some cases are then largely determined by their values in the other cases of the database where they are observable.

Integration of External Parameters in the Prior Distribution for Improved Estimation of Unobserved Parameters.

The key idea behind the use of priors is to model the distribution of parameters in the population. Then with the MAP estimation, the goal is to find the most likely parameters according to the correlations in this distribution. Here we explore the possibility of integrating parameters which are not estimated into the prior distribution (we call them external parameters), to influence the value of estimated parameters.

Namely, for all our 811 cases, the height and weight of the patients were available. We can also consider the heart rate which is not an estimated parameter. In order to add these three parameters to the prior distribution, we perform the following modifications to the method: in Equation 8, instead of considering a vector which only contains the parameters to be estimated, we use a concatenation vector of dimension 9 which contains the 6 estimated parameters and the 3 external parameters. Formulations of the prior covariance and mean in Equations 2 and in the equations of Section 2.3 are adapted as well to accomodate this concatenation vector.

We then perform the following simple experiment: from the estimated parameters in Section 4.3, we estimate the covariance matrix and the mean of the concatenation vector, which we then use as a prior for a new estimation. This is equivalent to perform only one iteration of the IUP algorithm with Full Matrix with this concatenation vector.

We report the most interesting result here in Figure 6 . For both the estimation without external parameters E1 (blue points) and with external parameters E2 (red points), we display the estimated values of the resting radius 0 for the cases where the Ejection Fraction is the only measurement, as a function of a/ the height of the patient (left), b/ the weight of the patient (middle) and c/ the Ejection Fraction of the patient.

The results are the following: first in both estimations the goodness of fit for the cases in the database is similar and high (a mean error of fit of around 0.06). For the estimation without external parameters E1, the resting radius is well correlated to the Ejection Fraction (this was already observed in Figure 5 ), but not at all to the weight and the height. However with the external parameters E2, the values of the resting radius is very correlated to the height and weight of the patients, and tends to increase with both these measurements. Since the resting radius, from a physical point of view, is related to the size of the heart, this correlation makes sense from a physical point

CONCLUSION

In this manuscript we introduced a method called Iteratively Updated Priors which performs successive personalisations of the cases in a population, where the prior probabilities on parameter values at each iteration are set from the statistics of personalised parameters in the previous iteration. We presented two different updates of prior parameters depending on the assumptions made on the distribution of parameters in the population, then derived a mathematical link between this algorithm and the minimisation of a population-wide cost function with a sparse penalty on the dimension of the set of personalised parameters. Depending on the type of covariance matrix update, this leads either to the selection of parameters which are set to a constant value in the population, or parameter directions which are set constant.

In Section 3 we demonstrated that the algorithm leads to a trade-off between the mean goodness of fit of personalised simulations over the database, and the dimension of the set of personalised parameters, and showed that this can be used to select a consistent subspace of minimal and sufficient dimension in which for each case of the population, there is a unique set of parameter values for which the simulation fits the corresponding measurements.

In Section 4, we presented an extended version of this algorithm applied to the personalisation of a larger database of 811 cases with heterogeneous measurements. This algorithms minimises the number of dimensions of the set of personalised parameters while ensuring a high (almost exact) goodness of fit over the whole database. We showed that the resulting population-based priors computed through this algorithm simulataneously leads to consistent parameter estimation for the cases where the parameters are not all observable, and are in particular influenced by their values in cases where they are observable. Finally, we show it is possible to add external parameters such as the height and the weight of the patients in the prior probability distribution in order to further guide the estimation of a consistent value of parameters which are unobservables in some cases.

A first extension is that beside the two types of updates presented here (Full Matrix and Diagonal Matrix), many other assumptions on the probability distribution of parameters in the population could be used, resulting in different formulations of the prior parameters. For example, the assumption could be made that the contractility 0 is independent of the other parameters, resulting in a block matrix covariance. Similarly, if the mean of a parameter value is known, then it is also easy to change the updates to accommodate this case (the derivations of the maximum likelihood updates for parameters of Gaussian priors are easily tractable even with such constraints). Finally, other type of prior shapes (other convex functions, mixture of gaussians...) could be used and the IUP algorithm would lead a dimensionality reduction of the parameter space and the creation of population-based statistics along the degrees of the freedom of the prior function (in our example it was the along the principal directions of the covariance matrix of the Gaussian prior).

Secondly, a last extension of this work would be to apply the IUP to a 3D cardiac model. However in the case of our 3D model [START_REF] Marchesseau | Fast parameter calibration of a cardiac electromechanical model from medical images based on the unscented transform[END_REF] whose simulations can be up to 1 hour of computation, the computational burden associated with performing a dozen or more repeated personalisations of the whole database can be too high. To tackle this and the 3D model, we believe a possibility would be to use the personalisation in a reduced subspace of 0D model parameters (such as H * ) to influence the estimation of a 3D model. Another possibility could be to investigate the use of specific 0D / 3D multifidelity parameter couplings as presented in [START_REF] Molléro | Multifidelity-cma: a multifidelity approach for efficient personalisation of 3d cardiac electromechanical models[END_REF] for multiple iterations of the IUP algorithm at once to lower the burden of repeated personalisations. We believe this could enable the use of the IUP algorithm on databases of 100 to 200 cases with the 3D model in around 4-5 days.

We conclude with a general comment on priors and the relevance of parameter selection in the context of personalisation. In the ideal case where we have prior statistics on all the parameters, every parameter which has an influence on a specific measurement will likely vary from the mean during parameter estimation, in accordance to the probability defined by the prior. This is why all the parameters vary at least slightly from the mean when only the ejection fraction or the stroke volume are available in our example. However if no statistics or external information is available on a set of parameters, a specific set of values of these parameters does not make more sense than another one that also fits the measurements. There is thus no other choice than arbitrarily choosing a direction, until more information can be given on this parameter direction from more experimental or physiological knowledge. Within these considerations, the best approach consists in our opinion in starting with a high number of parameters, and iteratively reducing the parameters space until a sufficient subspace is found in which all directions can be observed. We argue that our Iteratively Update Prior method, though it absolutely requires supervision for the selection of parameters with physiological relevance, gives a complete framework for consistent parameter estimation through the joint selection of a parameter subspace of reduced dimension and the creation of relevant priors on its directions.

Compliance with Ethical Standards : Informed consent was obtained from the subjects and the protocol was approved by the local Research Ethics Committee.

APPENDIX A: LINK BETWEEN THE IUP ALGORITHM AND THE ITERATIVELY REWEIGHTED LEAST SQUARE (IRLS) ALGORITHM

Here we demonstrate a link between the IUP algorithm and the Iteratively Reweighed Least Square in sparse regression. In particular, we exhibit a populationwide cost function with sparse regularisation. The sparse penalty correspond to the rank of the covariance matrix of the estimated set of personalised parameters, which has two different interpretation depending on the covariance matrix update.

For a given iteration of the IUP algorithm, we note X = ( 1 , ..., ) the concatenation of some vectors of parameter values , = 1.. for all the cases, and Ô = ( ̂ 1 , ..., ̂ ). We define the population-wide cost function S (X , Ô) as the mean of all the individual cost functions for all the = 1.. cases:

S (X , Ô) = 1 ∑ =1 ̂ ( , ̂ , , Δ ), S (X , Ô) = 1 ∑ =1 ( , ̂ ) + ∑ =1 ( ) .
The first sum of this expression is the sum of the data-fit terms over all the cases and is the same expression at each iteration. The second sum has a different expression at each step depending on the updates which are performed, which impact the behavior of the algorithm. We explicit here these expressions and explain the resulting behaviour of the IUP algorithm with these updates.

"Diagonal Matrix" updates.

At each iteration have:

∑ =1 ( ) = ∑ =1 ( -) (Δ ) -1 ( -),
where Δ is a diagonal matrix computed from the personalised parameters k-1 at iteration -1 with Equation 7. We can break down the sum along the = 1.. coordinates (corresponding to the N personalised parameters): we note Δ the j-th component of the diagonal of Δ which formula is:

Δ = 1 ∑ =1 (( k-1 ) -( k-1 ) ) 2 .
We then have:

∑ =1 ( ) = 1 ∑ =1 ∑ =1 (( ) -( ) ) 2 Δ , ∑ =1 ( ) = ∑ =1 1 ∑ =1 (( ) -( ) ) 2 Δ .
To understand the behaviour of our algorithm with these updates, we can look at the similarity between our updates and a classic method called Iteratively Reweighted Least Square, used to solve the minimisation of cost functions involving norms such as

* = || -( )|| = ∑ ( -( )) . ( 20
)
The IRLS algorithm works by succesively optimising the following (and usually easier) "weighted 2 " problems where in the k-th problem, each of the coordinates are reweighted by a vector derived from the norm of the coordinate in the solution of the previous problem:

( * ) = ∑ ( -( )) , where = | -(( * ) -1 )| -2 . ( 21
)
In this context, the Diagonal Matrix updates in our algorithm are reweighting the penalty of each coordinate (or personalised parameters) with the weight (Δ ) -1 which is inverse of the variance of the parameter in the current personalisation. To understand how the IRLS method applies in our case, we define the random vector as the difference ( -) of the personalised parameters to the mean of the population , for which all the cases are samples. The variance of each personalised parameter is then the ( 2 ) norm of the coordinate of this random vector. In this formulation, each update correspond to the IRLS formulation with p=0 for the minimisation of the 0 norm of .

As a consequence, if the successive and converge, according to the IRLS method the value of regularization term converges to:

∑ =1 ( * ) = || || 0 ,
where || || 0 is the number of coordinates of which are non-zero, which is also the number of personalised parameters which do not have a unique value in the population.

"Full Matrix" updates.

At each iteration we have:

∑ =1 ( ) = ∑ =1 ( -) (Δ ) -1 ( -).
As Δ is the covariance matrix of the -1 it can be orthogonally diagonalised and expressed as Δ = ( )( )( ) with an orthogonal matrix which columns are the principal directions of the set of -1 , diagonal which diagonal coefficients are the variances of the -1 in these principal directions. We can write:

∑ =1 ( ) = ∑ =1 ( ( -)) ( ) -1 ( -),
and in this case, each update consists in an IRLS step in the basis , defined by the principal directions of the set of -1 (the columns of ) in which the covariance matrix Δ is diagonal. Each of these updates thus correspond to an optimisation step of the 0 norm of the random vector of the difference ( -) of the personalised parameters to the mean of the population , expressed in the basis . As a consequence, if the successive and respectively converge to * and * , and we can build a random vector * on a basis * to which the and respectively converge, the regularization term converges to:

∑ =1 ( ) = || * || 0 ,
where || * || 0 is the number of coordinates of * which are non-zero expressed in the basis * (made of principal directions of the set of * ). It is the rank of X or more simply the number of principal directions in which parameters do not have a unique value in the population.

APPENDIX B: THE 0D CARDIAC MODEL

The 0D cardiac mdoel is a reduced version of our in-house 3D electromechanical model [START_REF] Marchesseau | Multiplicative jacobian energy decomposition method for fast porous visco-hyperelastic soft tissue model[END_REF][START_REF] Marchesseau | Fast parameter calibration of a cardiac electromechanical model from medical images based on the unscented transform[END_REF] , derived with the method presented by Caruel et al. [START_REF] Caruel | Dimensional reductions of a cardiac model for effective validation and calibration[END_REF] . In this section we first present the mechanical model itself, then we give the equations of the 0D model.

The Mechanical Model

The electromechanical equations are based on the Bestel-Clement-Sorine model (BCS) of sarcomere contraction for 3D models as extended by Chapelle et al. [START_REF] Chapelle | Energy-preserving muscle tissue model: formulation and compatible discretizations[END_REF] , in conjunction with a Mooney-Rivlin energy for the passive hyperelasticity. Hemodynamics are represented through global values of pressures and flows in the cardiac chambers, and coupled to the mechanical equations with the Windkessel model of blood pressure for the after-load (aortic pressure).

The BCS model: Active Contraction and Passive Material

The BCS model describes the sarcomere forces as the sum of an active contraction force in the direction of the fibre, in parallel with a passive isotropic visco-hyperelastic component. It is compatible with the laws of thermodynamics, and allows to model physiological phenomena at the sarcomere scale which translate at the macroscopic scale (such as the Starling Effect).

The active force in the sarcomere is modeled by the filament model of Huxley [START_REF] Huxley | Muscle structure and theories of contraction[END_REF] , which describes the binding/unbinding process of the actin and myosin in the sarcomere at the nanoscopic scale. At the mesoscopic scale, it results 3 in a differential equation which relates the active stress c , the stiffness c and the strain c of the filament within the sarcomere:

̇ c = -(| | + + | | -+ | ̇ c |) c + 0 | | + , ̇ c = -(| | + + | | -+ | ̇ c |) c + ̇ c c + 0 | | + , ( 22
)
where is a constant related to the cross-bridge destruction during contraction, 0 and 0 are respectively the maximum stiffness and contraction. The values of | | + and | | -are respectively the rate of build-up ATP and decrease RS of the force during contraction and relaxation, which depends on the depolarisation and repolarisation times and of the sarcomere:

= ⎧ ⎪ ⎨ ⎪ ⎩ ATP when ≤ ≤ -RS otherwise | | + = max( , 0), | | -= -min( , 0). (23)
This active force is applied in the direction of the fibre through the visco-elastic component, made of a spring and a dissipative term . As derived in Caruel et al. [START_REF] Caruel | Dimensional reductions of a cardiac model for effective validation and calibration[END_REF] , the resulting stress 1D in the fibre direction is given by:

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 1D = 1D -c (1 + 2 c ) 2 , ( + ̇ ) = ( 1D -c )(1 + 2 1D ) (1 + 2 c ) 3 , ( 24
)
where 1D = 1 ⋅ ⋅ 1 is the strain in the fibre direction 1 ( is the Green-Lagrange strain tensor).

Finally for the passive component the isotropic Mooney Rivlin model of hyperelastic material is used, driven by the following strain energy:

= 1 ( 1 -3) + 2 ( 2 -3) + 2 ( -1) 2 , ( 25
)
where 1 , 2 and are the invariants of the Cauchy-Green deformation tensor, 1 ,

2 and are the parameters of the material.

Haemodynamic Coupling

To model the influence of blood dynamics during the cardiac circle, the mechanical equations are coupled with a basic circulation model implementing the 4 phases of the cardiac cycle. For a given ventricle, if we note at the pressure in the atrium, ar the pressure in the artery and V the pressure in the ventricle, the phases are the following:

• Diastolic Filling: when V ≤ at , the atrial valve is open and the ventricle fills up with blood.

• Isovolumetric contraction: when contraction starts, V at ≤ V ≤ ar and all the valves are closed. • Systolic Ejection: when V ≥ ar , the arterial valve opens and the blood is ejected into the artery.

• Isovolumetric relaxation: when the contractile forces disappear, V finally decreases. at ≤ V ≤ ar again and all the valves are closed.

We use the haemodynamic model introduced by Chapelle et al. [START_REF] Chapelle | Energy-preserving muscle tissue model: formulation and compatible discretizations[END_REF] which links the blood flow to the ventricular, atrial and arterial pressures with the following equations:

= ⎧ ⎪ ⎨ ⎪ ⎩ at ( V -at )
for V ≤ at iso ( Vat ) for at ≤ V ≤ ar ar ( Vat ) + iso ( arat ) for V ≥ ar (26)

Here the atrial pressure at ( ) (cardiac preload) is imposed at a constant value at_lower except for a pressure bump up to at_upper at the beginning of cardiac cycle, to account for the contraction of the atrium before the ventricular contraction. Finally the pressure of the artery ar (cardiac afterload) is modeled with the 3-parameters Windkessel model [START_REF] Westerhof | Analog studies of the human systemic arterial tree[END_REF] and coupled to the ventricular outflow q through the equation:

̇ ar + ar -ve = ( + ) + ̇ , ( 27
)
where is the Peripheral resistance, is the Characteristic impedance, is the Arterial compliance and is the Central Venous Pressure.

The equations of the 0D model

We derived the equations of the 0D model of the heart from our 3D model with the method of Caruel et al. [START_REF] Caruel | Dimensional reductions of a cardiac model for effective validation and calibration[END_REF] . With this approach, both the 0D and the 3D model rely on the same BCS equations of myocardial forces, but the following simplifying assumptions are made on the geometry, the electrical activation and the properties of the material from the 0D model:

1. The ventricle has a spherical shape.

2. The material is incompressible.

3. The electrical activity is synchronous and homogeneous over the sphere.

With these assumptions, myocardial forces and motion are spherically symmetric and can be entirely described by the inner radius of the ventricle. Deformation and stress tensors can also be reduced to a simple form [START_REF] Caruel | Dimensional reductions of a cardiac model for effective validation and calibration[END_REF] , which leads to a system of a dozen equations which we report in Table 11 and Table 12 and detail in Section 7.2.1 and 7.2.2.

The equations are implemented into C code and the system of equations is solved using an explicit Forward Euler method with a temporal discretisation of 0.01 milliseconds. This leads to the simulation of 15 beats per second. Finally, the 0D model was also encoded in the CellML format [START_REF] Cuellar | An overview of cellml 1.1, a biological model description language[END_REF] , an open standard based on the XML markup language to store and exchange computer-based mathematical models. It can be downloaded from the Physiome Model Repository3 and easily exploited through the software OpenCOR 10 .

Mechanical Equations

The list of simplified equations of our 0D model is reported in Table 11 . Equations (a), (b), (c) and (f) are the same sarcomere and visco-elastic equations than Equations 22 & 24, which are calculated once for the whole sphere. in equations (d), (e), (g) and (h) denotes a component of the simplified Cauchy-Green deformation tensor which depends only on = -0 . passive in equation (g) is the stress due to the passive law and viscosity in equation (h) is the stress due to an additional viscous damping , both expressed as a simple function of (see [START_REF] Caruel | Dimensional reductions of a cardiac model for effective validation and calibration[END_REF] for the full derivations). In equation (i), Σ sph is the sum of all the stresses applied to the sphere. Equation (j) is the resulting equation of motion which, coupled with the haemodynamic model (k) and the windkessel equation (l), gives the full system of 3 equations to be solved at each iteration.

Electrophysiology Equations

Assuming synchronous and homogeneous electrical activation (and thus sarcomere force) means that all of the ventricle is depolarised simultaneously. This leads to a rate of ventricular pressure rise during the isovolumetric contraction (resp. isovolumetric relaxation) which is very close to the rate of build-up ATP (resp. decrease RS ) of the active stress . However in 3D, this rate is also very dependent on the time for the ventricle to be fully depolarised, which is roughly the QRS duration.

In order to correct this discrepancy between the models, we adapted the electrical parameter to take into account the QRS duration. We model the fraction depo of the ventricle which is currently depolarised as a piecewise linear function of time which depends on d,global , r,global and duration . Then the values of | | + and | | - in Equation (a) are adapted to depend on the value of depo as described in Table 12 .

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩

FIGURE 3

 3 FIGURE 3 Top: value of the eigenvalues of the SVD of personalised parameters at each IUP iteration. Bottom: mean values ̄ (X , Ô) and ̄ (X ) of the data-fit term ( , ̂ ) (blue) and the regularization term ( ) (green) across the population at each IUP iteration.

FIGURE 4

 4 FIGURE 4 Comparison of the values of the error of fit √ ( , ̂ ) in the FULL-IUP-05 and FULL-IUP * -05 runs.

FIGURE 5

 5 FIGURE 5 Top: (resp Bottom:) Log-value of estimated parameters value as a function of the Ejection Fraction (resp Stroke volume) for the 38 (resp 45) cases where only the Ejection Fraction (resp Stroke Volume) is available, in blue with the IUP* algorithm, in red without priors.

FIGURE 6

 6 FIGURE 6Estimated values of the resting radius 0 in the estimation without external parameters (blue points) and with external parameters (red points) in the prior distribution, for cases where only the Ejection Fraction was available (the resting radius was thus not completely observable).

  of view and leads to an improved consistency of the estimated values in E2 than in E1.

̇( 1 + 2 c ) 3 =) 2 passive = 4 ( 1 -

 13241 c = -(| | + + | | -+ | ̇ c |) c + 0 | | + ̇ c = -(| | + + | | -+ | ̇ c |) c + ̇ c c + 0 | | + ( + ̇ ) = ( 1Dc )(1 + 2 1D ) -3 )( 1 + 2 ) viscosity = 4 (1 + -6 ) ̇ Σ sph = 1D + passive + viscosity 0 ̈ = (1 ( Vat ) for V ≤ at iso ( Vat )for at ≤ V ≤ ar ar ( Vat ) + iso ( arat ) for V ≥ ar ̇ ar + arve = ( + ) + ̇

TABLE 1

 1 Sets of 0D model parameters and outputs for which the measurements are available in the example.

	Personalised Parameters P
	Contractility 0
	Resting Radius 0
	Stiffness 1
	Peripheral Resistance
	Windkessel Characteristic Time
	Venous Pressure ve
	Outputs O
	Maximal Volume
	Stroke Volume
	Mean Aortic Pressure
	Diastolic Aortic Pressure
	Name		Update
	FULL-IUP-005 0.05
	FULL-IUP-01	0.1	Full Matrix
	FULL-IUP-05	0.5
	DIAG-IUP-005 0.05
	DIAG-IUP-01	0.1	Diagonal Matrix
	DIAG-IUP-05	0.5

TABLE 2

 2 Names and parameters of the 6 runs of the IUP algorithm

TABLE 3

 3

	Run		≥ 10 -4	̄ (X , Ô)	̄ (X )
	FULL-IUP-05	0.5	3	1.85	1.5
	FULL-IUP-01	0.1	4	2.15e-2	0.40
	FULL-IUP-005 0.05	4	6.18e-3	0.20
	DIAG-IUP-05	0.5	4	0.65	2.04
	DIAG-IUP-01	0.1	4	0.23	0.41
	DIAG-IUP-005 0.05	5	1.96e-2	0.25
	TABLE

For each run of the IUP algorithm, final eigenvectors of the SVD of personalised parameters ( ) and their corresponding eigenvalue . In bold we emphasize the coordinates of eigenvectors which are lower than 10 -3 and the eigenvalues which are lower than 10 -4 .

TABLE 5

 5 Mean error of fit and variability of personalised parameters in 5 personalisations for the parameter subspaces. The value is reported as 0 if it is lower than 10 -3 .

			√
				0	0	1	ve
	H 0	0	0.03	0	0.13	0.02	0.03 0.04
	H	1	0	0	0	-	0.02	0.02 0.03
	H ( 1 , 0 ) 0.99	-	0	-	0.02	0.02 0.14
	H ( 1 , 0 ) 3.90	0	-	-	0.05	0.02 0.04
	H ( 1 , ve ) 0.06	0	0	-	-	0	0.06
	H ( 1 , ) 0.07	0	0	-	0.07	-	0
	H ( 1 , )	0.01	0	0	-	0.02	0	-
	H 5	0	0	0	0	0.015 0.017 0.03
	H 4	0	0	0	0	0	0	0
	H 3	1.99 0.03	0	0.002	0	0.004 0.01
	H 2	3.03	0	0	0	0	0	0

TABLE 6

 6 Final eigenvectors and eigenvalues of the IUP run with Full Matrix and = 0.05 on H 5 .

			Final-IUP-005
		0	0	ve
	1	0.58 0.07	0.05 -0.51 0.63 4.77467
	2	0.46 -0.32 0.50	0.65	0.10 3.41794
	3	0.65 -0.06 -0.44 -0.06 -0.61 2.95973
	4	0.16 0.94	0.15	0.25 -0.05 1.03746
	5	0.01 -0.01 0.73 -0.50 -0.47 4.30e-6

TABLE 7

 7 Mean error of fit and variability of personalised parameters in 5 personalisations in H * . The value is reported as 0 if it is lower than 10 -3 .

	√
	0	0	ve

TABLE 8

 8 To build this database, we focused Measurements considered in the heterogeneous database on gathering patients and acquisitions for which the heart rate and at least one of the 6 following measurements in Table8was available:Within this database of 811 cases, we have the following statistics:

	Measurements O
	Maximal Volume
	Minimal Volume
	Stroke Volume
	Ejection Fraction
	Mean Aortic Pressure
	Diastolic Aortic Pressure

• Pressure Measurements are available for 651 cases only. • The Maximal Volume and Minimal Volume are both available for only 340 cases. • Among the 471 other cases, either the Ejection Fraction (63 cases), Stroke Volume (386 cases) or no volume measurement at all (21 cases) are available. • Ejection Fraction is the only measurement available in 38 cases. • Stroke Volume is the only measurement available in 45 cases. • 258 cases have a 'complete' set of 4 measurements: Maximal Volume, Minimal Volume, Mean Aortic Pressure and Mean Diastolic Pressure (we do not report Ejection Fraction and Stroke Volume if the Maximal Volume and Minimal Volume are already reported).

  and a corresponding normalisation vector N . The normalisation coefficient for the Ejection Fraction (EF) (which is a percentage) is 5%, it is 10 ml for the volume values (MaxV,MinV,SV) and 200Pa for the pressure values (MP, SP and DP).

	1. (	,		,	,	),	7. (	,	),
	2. (	,	,	,	),		8. (	,	),
	3. (	,	,	,	),		9. (	,	),
	4. (	,	,	),		10. (	),
	5. (	,	,	),		11. (	),
	6. (	,	,	),		12. (	),

TABLE 10

 10 Percentage of cases in the population with the corresponding final value of in Equation11at convergence of the FULL-IUP * -05 run.

	1	2	4	8	16
	% 81.79 0.49 1.47 9.85 6.4

TABLE 11

 11 Mechanical equations of the 0D model

	=	⎧ ⎪ ⎪ ⎪ 1 ⎨ ⎪ (1 --, duration -,	when , when , ) when ,	≤ ≤ , + duration ≤ ≤ , + duration ≤ ≤ , +
		⎪	duration
		⎪
		⎩

duration 0 otherwise | | + = ATP . depo | | -= RS .(1depo )

(28)

TABLE 12

 12 Electrical activation in the 0D model

A practical derivation of the maximum likelihood estimation of the covariance matrix can be found in https://people.eecs.berkeley.edu/ jordan/courses/260-spring10/other-readings/chapter13.pdf

https://models.physiomeproject.org/e/470

Funding:

The research leading to these results has received European funding from the FP7 project MD-Paedigree (600932), the ERC starting grant ECSTATIC (715093) and the ERC advanced grant MedYMA (2011-291080) and was supported by the grant AAP Santé 06 2017-260 DGA-DSH and by the Inria Sophia Antipolis -Méditerranée NEF computation cluster.

Conflict of Interest:

The authors declare that they have no conflict of interest.