# Inverse Uncertainty Quantification for Personalized Biomechanical Modeling: Application to Pulmonary Poromechanical Digital Twins

A. Peyraut, M. Genet

► To cite this version:

A. Peyraut, M. Genet. Inverse Uncertainty Quantification for Personalized Biomechanical Modeling: Application to Pulmonary Poromechanical Digital Twins. Journal of Biomechanical Engineering, 2025, 147 (8), pp.081003. 10.1115/1.4068578. hal-05118573

HAL Id: hal-05118573

https://hal.science/hal-05118573v1

Submitted on 18 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Inverse uncertainty quantification for personalized biomechanical modeling: application to pulmonary poromechanical digital twins

A. Peyraut¹,² and M. Genet¹,²

¹Solid Mechanics Laboratory, École Polytechnique/IPP/CNRS, France

²MΞDISIM Team, INRIA, France

June 17, 2025

## 1 Abstract

The development of personalized models is a key step for addressing various problems, especially in biomechanics. These models typically include many constants, introduced in the model material law or loading definition, and their estimation is crucial for the model personalization. However, performing solely the estimation does not yield any information on the estimation accuracy. Additionally, all parameters can typically not be estimated based only on clinical data: some parameters are identified, while others are fixed at generic values. The question of the identifiability of the parameters, along with the robustness of the estimation, notably to measurement errors and to model errors, is therefore crucial and should be quantitatively addressed in parallel to the model development. In this paper, we propose a general inverse uncertainty quantification pipeline based on the creation of synthetic data –for which the parameters ground-truth values are known–, generated for different noise and model error levels. Estimation is then performed for many realizations of the noise or model errors, as well as parameter initializations, until convergence of the estimated parameters error distributions. This pipeline was applied to a poromechanical lung model for illustration and validation purposes. It provides quantitative information on the actual identifiability of the parameters, and any derived quantity of interest, in the clinical setting. In particular, it allows to retrieve a confidence interval for each estimated parameters, which represents valuable information for diagnosis or prognosis use of the estimated values. This work is therefore a step towards improving the reliability of digital twins pipelines.

## 2 Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>Methods</td><td>4</td></tr><tr><td>2.1</td><td>Parameter estimation problem and inverse uncertainty quantification approach</td><td>4</td></tr><tr><td>2.1.1</td><td>General parameter estimation problem formulation</td><td>4</td></tr><tr><td>2.1.2</td><td>Proposed inverse uncertainty quantification approach</td><td>5</td></tr><tr><td>2.2</td><td>Brief summary of the lung poromechanical model</td><td>7</td></tr><tr><td>2.2.1</td><td>Input data</td><td>7</td></tr><tr><td>2.2.2</td><td>Loading and unloading problems</td><td>9</td></tr><tr><td>2.3</td><td>Brief summary of the model personalization pipeline</td><td>9</td></tr><tr><td>2.4</td><td>Implementation details</td><td>11</td></tr><tr><td>3</td><td>Results</td><td>13</td></tr></table>
<table><tbody><tr><td>3.1</td><td>Identification errors with one zone</td><td>13</td></tr><tr><td>3.1.1</td><td>Impact of measurement errors on the identified stiffness</td><td>13</td></tr><tr><td>3.1.2</td><td>Impact of measurement errors on the identified pleural pressure</td><td>13</td></tr><tr><td>3.1.3</td><td>Simultaneous estimation of the stiffness and of the pleural pressure</td><td>14</td></tr><tr><td>3.2</td><td>Identification errors with two zones</td><td>16</td></tr><tr><td>3.2.1</td><td>Analysis of [Patte et al. 2022a]'s parametrization</td><td>16</td></tr><tr><td>3.2.2</td><td>Analysis of [Laville et al. 2023]'s parametrization</td><td>16</td></tr><tr><td>3.3</td><td>Estimation with n zones</td><td>19</td></tr><tr><td>3.3.1</td><td>Estimation for a lung with 3 zones</td><td>19</td></tr><tr><td>3.3.2</td><td>Estimation for a lung with 4 zones</td><td>20</td></tr><tr><td>3.4</td><td>Improving the estimation process using gravity</td><td>20</td></tr><tr><td>3.4.1</td><td>Estimation of all stiffnesses and of the pleural pressure</td><td>21</td></tr><tr><td>3.4.2</td><td>Estimation with n zones</td><td>21</td></tr><tr><td>4</td><td>Discussion</td><td>22</td></tr><tr><td>4.1</td><td>Identifiability and sensitivity analysis</td><td>23</td></tr><tr><td>4.1.1</td><td>Parameter identifiability</td><td>23</td></tr><tr><td>4.1.2</td><td>Sensitivity analysis</td><td>24</td></tr><tr><td>4.1.3</td><td>Analysis of former parametrizations</td><td>25</td></tr><tr><td>4.2</td><td>Improving the estimation process by enriching the input dataset</td><td>26</td></tr><tr><td>5</td><td>Conclusion</td><td>26</td></tr><tr><td>References</td><td></td><td>33</td></tr><tr><td>Appendix</td><td></td><td>34</td></tr><tr><td>A</td><td>Recall of the equations of the model</td><td>34</td></tr><tr><td>A.1</td><td>Poromechanical framework</td><td>34</td></tr><tr><td>A.2</td><td>Equilibrium equations</td><td>35</td></tr><tr><td>A.3</td><td>Constitutive framework</td><td>37</td></tr><tr><td>B</td><td>Using gravity to improve the estimation process: a proof of concept on the linear elastic spring</td><td>39</td></tr><tr><td>C</td><td>Recall of the weak formulation of the problem</td><td>42</td></tr></tbody></table>

# 1 Introduction

In recent years, *in-silico* models have gained significant attention in the field of biomechanics [Hu et al. 2001; Fetita et al. 2019; Gonsard et al. 2024]. These computational models could indeed represent valuable tools for improving diagnosis, classification, prognosis, and patient care [Smith et al. 2011; Genet et al. 2014; Kaul 2019; Carniel et al. 2020]. Specifically, they could help better understand specific diseases [Aguado-Sierra et al. 2011; Genet et al. 2016; Agustí et al. 2022; Pellegrino et al. 2023], explore drug delivery mechanisms, and even assist surgeons in preoperative planning and practice [Lee et al. 2014; Rausch et al. 2017].

Given their clinical applications, these models must have a high degree of accuracy in simulating physiological behaviors. However, physiological systems vary significantly across patients, with factors such as anatomy and mechanical responses differing based on age, gender or medical history, for example [Hofmann et al. 2002; Sicard et al. 2018]. To develop personalized models, it is hence essential to incorporate patient-specific clinical data. Clinical data indeed provides critical information, such as anatomical geometry, which can be directly integrated into the model. Furthermore, this data also plays a crucial role in parameter estimation, a fundamental step in building patient-dependent models.
Parameter estimation is indeed one of the main steps of model personalization. Biomechanical models usually include many constants, introduced, for example in the formulation of the mechanical laws, in the loading, or in the parameterization of the system shape. If some parameters are likely to be similar across patients –e.g., healthy lung compliance [Mraz et al. 2024]–, other parameters may vary greatly [Cook et al. 2014; Plantier et al. 2018]. In order to describe accurately the behavior of a system, strongly patient-dependent parameters should therefore be estimated based on clinical data. This identification process, which corresponds to an inverse problem [Tarantola 2005], consists in estimating model parameters with computer methods and in studying the influence of input parameters on the global response of a model [Le Gratiet et al. 2014]. It offers many advantages, such as accessing to values that cannot be experimentally assessed, or only obtained with invasive procedures. For example, for pulmonary applications, parameter estimation could allow us access to pleural pressure –a subatmospheric pressure keeping lungs inflated– values, which are currently only accessible using esophageal balloons, which represents a very invasive procedures [Milic-Emili et al. 1964].

However, for complex models, all parameters of interest can usually not be estimated based solely on clinical data. Parameter non-identifiability can be inherent to the model, which can be assessed through “a priori” (or “structural”) identifiability approaches [Siekmann et al. 2012; Anstett-Collin et al. 2020; Wieland et al. 2021], or can be due to insufficient data, which can be studied with “practical” identifiability analyses [Daly et al. 2018]. Some parameters therefore need to be chosen at reference values, which are selected within a physiological range, and are hence potentially wrong, as variations around the physiological values can be expected across individuals. This necessarily induces an error for parameter estimation. Additionally, clinical data are typically noisy, which also introduces a bias in the estimation process. The chosen estimation method can also induce even more uncertainty in the estimation: many methods exist, which all display different computational costs and accuracies [Avril et al. 2008; Roux et al. 2020]. Assessing the structural identifiability of the parameters of the model, as well as investigating biases on the estimation introduced by noise and model errors, i.e., the practical identifiability, is therefore crucial in many biomechanical problems.

It is, however, impossible to know the exact error committed on the identified parameters during estimation processes, as their ground-truth values are unknown. To address this issue, various uncertainty quantification approaches have been proposed to evaluate the estimation robustness. If “direct” uncertainty quantification approaches study model-related uncertainties, typically involving sensitivity analyses, such as the Sobol approach [Saltelli et al. 1995; Sobol et al. 2007; Wirthl et al. 2023], to quantify the sensitivity of the model output to different input parameters, “inverse” uncertainty quantification approaches focus on evaluating uncertainties related to model parameter estimation [Wu et al. 2018]. These approaches can include Bayesian approaches [Ying 1991; Zhang 2004; Nagel et al. 2016; Kerfriden et al. 2019], which provide a probabilistic error distribution for the estimated parameters. The distribution is computed by assuming the form of the prior, i.e., by providing a probabilistic law describing the input parameters. These uncertainty quantification methods may, however, be too expensive for complex problems. In this article, we propose a statistical inverse uncertainty quantification pipeline applicable to computationally expensive problems that we validate on the *in-silico* lung model proposed by [Patte et al. 2022b], with the boundary conditions proposed in [Peyraut et al. 2024a].

[Patte et al. 2022b] developed a lung model to investigate Idiopathic Pulmonary Fibrosis (IPF), a specific Interstitial Lung Disease (ILD) [Antoniou et al. 2014]. ILDs affect the mechanical properties of the lung parenchyma, in particular through abnormal collagen deposits [Bates et al. 2007; Thannickal et al. 2014; Wijsenbeek et al. 2020]. Investigating ILDs, and hence IPF with a mechanical framework is therefore particularly interesting. IPF is an idiopathic disease; the mechanisms driving its appearance and evolution are therefore poorly understood. However, a hypothesis suggesting the existence of a mechanical vicious cycle driving the progression of IPF [Hinz et al. 2016; Wu et al. 2020] has been proposed, making the mechanical study of the disease particularly interesting. According to this hypothesis, abnormal collagen deposits occurring in IPF during tissue repair processes, which leads to stiffer lung tissues, cause higher mechanical
stresses during the breathing cycle. These stresses being non-physiological, fibroblasts, the cells responsible for collagen production, are stimulated [Murata et al. 2014; Hall et al. 2024], which in turn leads to further collagen deposition and fibrosis, and hence higher tissue stiffness.

As a result, a digital lung model for patients with IPF could be particularly useful. It could in particular help investigate correlation of mechanical biomarkers, such as stress distributions or reduced lung compliance, with fibrosis progression. Such a model could help clinicians with prognosis and improve patient care. The development of a digital twin relies heavily on the available clinical data [Bordas et al. 2015; Roth et al. 2017; Coorey et al. 2022]. In the case of IPF, French clinical routine imaging typically includes two CT scans –one at end-exhalation and one at end-inhalation–, from which a porosity field can be extracted. To integrate this clinical data, and hence improve the personalization, [Patte et al. 2022b; Peyraut et al. 2024a] developed a quasi-static poromechanical model, recalled in this article.

In this paper, we describe a statistical inverse uncertainty quantification pipeline, and apply it to the developed pulmonary model. We use this pipeline to perform a quantitative study of the identifiability (both structural and practical) of the parameters of our model, by studying the estimation errors for different parameterizations, i.e., for different sets of parameters identified. As the ground-truth values of the parameters is unknown for clinical applications –i.e., for a real patient–, we performed our study using synthetic data. Additionally, as the identification might be very sensitive to noise and/or to model errors [Avril et al. 2004; Rossi et al. 2015], we investigate the impact on the estimation when such errors are added to the generated synthetic data.

This paper is organized as follows. In the Methods Section, we present our general inverse uncertainty quantification pipeline, and present the application to [Patte et al. 2022b]’s poromechanical framework. Additionally, the impact of the input data on the estimation quality is discussed. Next, in the Results section, we propose the analysis of different parameters parameterizations, by studying their noise and model errors robustness. The best parameterization to choose based on the identified model sensitivity and based on the noise and model errors robustness is then thoroughly discussed in the Discussion Section.

# 2 Methods

## 2.1 Parameter estimation problem and inverse uncertainty quantification approach

### 2.1.1 General parameter estimation problem formulation

Estimation of model parameters, e.g., material and loading constants, from experimental data is a central step towards model personalization. Many different approaches have been proposed to tackle inverse problems when full-field measurements are available. The most common are the FEMU method [Kavanagh et al. 1971], the Equilibrium Gap Method [Claire et al. 2004] and the Virtual Field Method [Grédiac et al. 2006]. An overview of the different methods and their performances can be found in [Avril et al. 2008] and [Roux et al. 2020].

All these methods rely on building a cost function, denoted by $L$ in this paper, which takes as input the measured displacement field, which we call $\underline{U}_{meas}$, and depends on the parameters to identify. The estimated parameters are denoted as $\underline{\theta}$ in this work, while the parameters fixed at reference values during the estimation are denoted as $\underline{\eta}$. The estimation then simply consists in finding the optimal parameters $\underline{\theta}$, denoted as $\underline{\theta}^{opt}$, which minimize the cost function $L$. The problem can therefore be written

$$ \underline{\theta}^{\text{opt}} := \underset{\underline{\theta}}{\operatorname{argmin}} \{L(\underline{U}_{\text{meas}}, \underline{\theta}, \underline{\eta})\}. \qquad (1) $$

For complex problems, it is however rare to know if the set of chosen parameters $\underline{\theta}$ is identifiable without prior studies. In the case of biomedical applications, estimating quantitatively the reliability of our identified parameters is however a key issue, since biomarkers could be derived from these
quantities. It is therefore necessary to assess the relevance of the chosen parametrization -i.e., the choice of the parameters to identify $\underline{\theta}$ and the choice of the parameters to fix at reference values $\underline{\eta}$- and to have precise information on the error committed during the estimation, especially as the clinical data used as input is often noisy, and as fixing the parameters $\underline{\eta}$ at reference values may introduce an error on the estimation.

We therefore propose in this article the formulation of an inverse uncertainty quantification pipeline, allowing us to study the confidence level of the estimation of the chosen parameters. This pipeline notably allows us to identify relevant parametrizations, along with studying the robustness of the estimation to measurement and model errors. For illustration purposes, our pipeline is then applied to pulmonary poromechanics but could easily generalize to any inverse problem.

### 2.1.2 Proposed inverse uncertainty quantification approach

In order to quantify uncertainty on estimated parameters, noisy and/or biased synthetic data should be created. The parameters $\underline{\theta}$ are then estimated using Equation (1) with the created synthetic data. If the estimation was performed with actual clinical data, it would indeed not give information on the error committed during the estimation.

**Synthetic data** To generate synthetic data, we begin with choosing the identified parameters $\underline{\theta}$ and the fixed parameters $\underline{\eta}$ at physiological values $\underline{\theta}_0$ and $\underline{\eta}_0$ -the ground-truth value of $\underline{\theta}$, $\underline{\theta}_0$, is hence known-. The associated displacement field $\underline{U}_0 = \underline{U}(\underline{\theta}_0, \underline{\eta}_0)$ is then computed using a model. Clinical images typically include noise, which has an impact on the estimation quality. To take this phenomenon into account, a random Gaussian noise, $\underline{f}$, with null mean and a standard deviation fixed according to a given level of noise, is added to the synthetic displacement field. The level of noise is defined by the Signal-to-Noise Ratio, a.k.a. SNR. As a result, our synthetic measure, $\underline{U}_{noisy}$, is defined as:

$$ \underline{U}_{noisy} := \underline{U}_0 + \underline{f}. \qquad (2) $$

As mentioned previously, it is rarely possible to identify all parameters of a model. As a result, during the estimation, the parameters $\underline{\eta}$ are fixed at reference values, chosen within a physiological range. However, this necessarily induces a bias in our model, since the values chosen for the reference parameters do very likely not correspond to the ground-truth values. These reference values are indeed chosen within a physiological range, while those values may vary greatly from a patient to another; it is for example the case for the lungs properties [Patte et al. 2022a].

The error introduced by the fixed parameters, which we refer to as “model errors” and corresponds to epistemic errors, can be modeled by performing the estimation with “wrong” values of $\underline{\eta}$. This can be achieved by adding an error $\underline{b}$ to the ground-truth value $\underline{\eta}_0$. As a result, the fixed parameters used during the estimation $\underline{\eta}_{epistemic}$ are defined as:

$$ \underline{\eta}_{epistemic} := \underline{\eta}_0 + \underline{b}. \qquad (3) $$

**Estimation** The optimized parameters are then estimated by optimizing the cost function presented in Equation (1). The estimated parameters are initialized randomly in a normal law centered on the ground-truth value and with a standard deviation chosen arbitrarily at 30 % of the ground-truth value. This random initialization allows us to account for local minima, which can bias the estimation: if the output distribution is the same as the input distribution, numerous local minima exist and the estimation is not possible, which would not be detected if the parameters were initialized at their ground-truth value. The estimation pipeline is described in Figure 1.

**Error quantification** The noise added to the synthetic displacement field is generated randomly. As a result, for a given noise level, the estimation process is repeated for many realizations of the noise, until convergence of the error distribution on the identified parameters is reached for the said noise level. The convergence of the error distribution is defined as follows.
<figure><img src="image_2.png" /><figcaption>Figure 1: Description of the proposed inverse uncertainty quantification pipeline. Synthetic displacement fields $U_{\text{noisy}}$ are generated by adding noise to a reference displacement field $U_0$, computed with the model, when the model parameters are chosen at reference values. Model errors can also be introduced by adding noise to the fixed parameters $\eta$. The cost function is then minimized using the derivative-free algorithm Nelder-Mead.</figcaption></figure>

Let us consider N converged processes. If the error distribution including N converged values of the estimated parameters $\theta$ is denoted $D_{\theta}(N)$, the distribution is converged if the $N+1$ value verifies

$$ \max_{\mathcal{F} \in [\sigma, \mu, q_1, q_2]} \left| \frac{\mathcal{F}(D_{\theta}(N+1)) - \mathcal{F}(D_{\theta}(N))}{\mathcal{F}(D_{\theta}(N))} \right| < \text{tol}, \quad (4) $$

where tol is the tolerance, chosen arbitrarily, for considering that the distribution converged. $\mathcal{F}$ represents different measures: the standard deviation $\sigma$, the mean of the distribution $\mu$, the first quartile $q_1$ or the third quartile $q_3$. In practice, we impose N to have a minimal value, i.e., we only check the convergence after a minimum number of iterations. This convergence criterion therefore ensures that we have a certain statistical stability on our error distribution for it to be converged: any new added value should not change significantly the distribution.

**Identifiability** In this work, we study the identifiability of model parameters. As discussed in [Siekmann et al. 2012; Anstett-Collin et al. 2020; Wieland et al. 2021], we can distinguish “a priori” or “structural” identifiability from “practical” identifiability (which considers measurement). If we do not conduct a theoretical identifiability analysis, our approach allows us, however, to study both types of identifiability computationally and quantitatively. We indeed assess parameter identifiability by first generating a reference parameter-to-data map using the model, for reference
parameters values. We then attempt to invert this map by using an optimization process. We do this many times, each time starting with a random initial value in a neighborhood of the solution. In case the final distribution of the estimation error (i.e., the distance between estimated and reference values) tends (in a statistical sense) towards a Dirac distribution centered on zero, we then claim the parameters are structurally identifiable, since the monotony of the parameter-to-data map is empirically established. We then add noise to the reference parameter-to-data map, and repeat the optimization process. Again, we run the estimation many times, for many realizations of the noise and many initializations of the optimization. This process assesses the practical identifiability of the parameters. In this case, if the final distribution of the estimation error tends towards a Dirac distribution centered on zero in the absence of noise, and if the error committed on the estimated parameters is reasonably small for reasonable noise levels, we claim our parameters are practically identifiable [Daly et al. 2018].

## 2.2 Brief summary of the lung poromechanical model

We illustrate our proposed inverse uncertainty quantification pipeline on our recently proposed lung digital twin pipeline [Patte et al. 2022a; Patte et al. 2022b; Laville et al. 2023]. The model depends on clinically available data, detailed hereafter.

### 2.2.1 Input data

In this article, we investigate the identifiability of different parameterizations, supposing that we have access to full-field experimental measurements. Note that the experimental measurements can differ depending on the investigated problem. The inverse uncertainty quantification pipeline presented here, however, requires at least a displacement field between at least two images. The information may be richer, e.g., more than one experimental measurement is available, which typically leads to a better estimation quality.

**Standard clinical images** The clinical dataset we have access to in the scope of our pulmonary digital twin consists in two CT-scans in supine position, one at end-exhalation and one at end-inhalation. This case corresponds to routine imaging practice in France for pulmonary fibrosis applications. In order to be patient-dependent, we base our model on the available data. As only two images are available, we choose to build a quasi-static lung model. Additionally, in order to include as much data as possible, and as we can extract a porosity field from the images, we choose to work in a poromechanical framework. We then propose to study the estimation quality of the parameters of this model, using in a first approach solely images in supine positions, as it corresponds to our real clinical problem.

**Additional images** The image dataset could however be enriched by considering additional images, e.g., images in prone position. In this case, the input data would be composed of up to four displacement fields: the displacement field from end-exhalation to end-inhalation in prone and in supine position, and the displacement field from prone to supine position at end-exhalation and at end-inhalation. Displacement fields from images in prone and supine positions are used because the end-exhalation to end-inhalation displacement differs depending on the patient's position [Barrau 2024]. In a first approximation, we can consider that the only difference between the two positions is a change in gravity orientation, which is a perfectly known force included in our model. As a result, we also investigate in this article whether including images in both supine and prone positions can improve the estimation. The principle is described in Figure 2. A proof of concept on a very simple case –a linear spring–, highlighting the potential advantages of enriching the dataset by using gravity, is presented in Appendix B.

Additionally, if CT is the traditional imaging modality for pulmonary applications, images can also be obtained with MRI techniques. In this case, the data set is very rich, and the input is composed of many displacement fields, the breathing cycle being composed of many intermediate images. An example of different possible datasets for modeling a breathing cycle is presented in Figure 2.
<figure><img src="image_3.png" /><figcaption>Figure 2: Illustration of different datasets for the estimation with a model representing a breathing cycle. Only the prone and supine positions are considered. $\underline{U}$ represents the displacement field between two positions, the superscripts p and s refer to the prone and the supine position respectively. Additional information can be obtained by adding new displacement fields as input, e.g., by adding CT-scans in another position. Additionally, the subscript ti (i = 0, 1, ..., n) refers to an image at a time ti, when n MRI images are considered. The estimation can be, for example, performed with one displacement field as input, which is usually the case with CT-scans. When MRI images are available, a large set of input data is available.</figcaption></figure>

In general, the estimation quality depends largely on the available dataset. We can expect a better estimation when more data is available. More images indeed contain more information, and can help to mitigate the effects of measurement errors, for example.

**Motion tracking** In order to extract the displacement field between two images, motion tracking is usually used [Grédiac et al. 2013; Sutton et al. 2015]. Motion tracking consists in finding the displacement field between an image $I_0$ and an image $I_t$ by minimization of an image energy $\Xi^{ima}$, to which a regularization energy $\Xi^{regul}$ is added to force a mechanically consistent solution [Leclerc et al. 2011; Mendoza et al. 2019]. The tracking problem can be written

$$ \underline{U}_{\text{meas}} := \underset{\underline{U}}{\operatorname{argmin}} \left\{ (1 - \beta) \Xi^{\text{ima}}(\underline{U}) + \beta \Xi^{\text{reg}}(\underline{U}) \right\}, \quad (5) $$

where $\beta$ represents the regularization strength.

The expression chosen for the image energy differs given the application [Škardová et al. 2019]. Lung specific image correlation processes have been developed in [Hurtado et al. 2017; Boucneau et al. 2020; Maghsoudi-Ganjeh et al. 2021]. In this article, we choose to use the following classical expression:

$$ \Xi^{\text{ima}}(\underline{U}) := \frac{1}{2} \int_{\Omega} (I_t(\underline{X} + \underline{U}(\underline{X})) - I_0(\underline{X}))^2 \, d\Omega. \quad (6) $$

Similarly, different expressions of the regularization energy have been proposed [Miller et al. 1993; Veress et al. 2005]. In this article, we will use the equilibrium gap regularization introduced in [Genet
## 2.2.2 Loading and unloading problems

**Problem description** The “measured” displacement we have access to consists in the displacement field between end-exhalation and end-inhalation. The computation of this displacement field with our model is, however, not straightforward: the end-exhalation configuration is indeed already loaded, in particular with the end-exhalation pleural pressure and gravity. The computation of the displacement field therefore requires two steps: an unknown unloaded configuration is first computed from the known loaded end-exhalation configuration. The unknown deformed end-inhalation configuration is then computed from the previously computed unloaded configuration. Note that by assuming that the unloaded configuration is stress-free, we implicitly neglect residual stresses, which could however be introduced when experimental evidence becomes available [Genet et al. 2015a; Genet 2019]. The problem is detailed in Figure 3. The poromechanical framework, equilibrium equations and associated notations are reminded in Appendix A. For illustrative purposes, Figure 4 presents the finite-element solution of the problem solved.

<figure><img src="image_4.png" /><figcaption>Figure 3: Steps of the computation of a breathing cycle with our poromechanical model: (a) unloading problem, i.e., computation of the unknown unloaded configuration from the known end-exhalation configuration, (b) loading problem, i.e., computation of the unknown end-inhalation configuration from the known unloaded configuration. All operators appearing in this Figure are formulated in Appendix A.</figcaption></figure>

## 2.3 Brief summary of the model personalization pipeline

As mentioned in the Introduction, building a personalized *in-silico* model can provide clinicians with a very interesting tool for understanding better some diseases and help them improve patient care.

**Estimation method** In the scope of the previously presented problem, the chosen estimation method is the Finite Element Model Updating (FEMU). Even though expensive -it requires the computation of the model at each iteration-, this method is extremely robust and displays a low noise and model error sensitivity. The FEMU method relies on the very intuitive principle of minimizing the gap between the measured displacement field and a displacement field computed with our model, which depends on the different parameters. The optimized parameters correspond
<figure><img src="image_5.png" /><figcaption>Figure 4: Illustration of the problem solved for with our poromechanical pulmonary model.</figcaption></figure>

$$L^{\text{FEMU}}(\theta, \eta) := \frac{1}{2} \|\underline{U}_{\text{meas}} - \underline{U}(\theta, \eta)\|_{L_2(\Omega)}^2 \\ \approx \frac{1}{2} \|\underline{U}_{\text{meas}} - \underline{U}(\theta, \eta)\|_{L_2}^2, \qquad (7)$$

$$\|f\|_{L_2} = \sqrt{\int |f(x)|^2 dx}, \qquad (8)$$

$$\|f\|_{l_2} = \sqrt{\sum_{k=1}^{n} |f_k|^2}, \qquad (9)$$

$$\frac{\theta_i - \theta_{i0}}{\theta_{i0}} \times 100. \qquad (10)$$

<footnote n="1">1 to the parameters allowing us to find a computed displacement field as close as possible to the measure.</footnote>

<footnote n="3">3 The cost function of the FEMU method can be formally expressed as</footnote>

<footnote n="4">4 with $\underline{U}_{\text{meas}}$ and $\underline{U}(\theta, \eta)$ is the nodal vector representing the finite element discretizations of $\underline{U}_{\text{meas}}$ and of $\underline{U}(\theta, \eta)$ respectively. In this case, the norm $L_2$ represents</footnote>

<footnote n="6">6 while $l_2$ represents</footnote>

<footnote n="7">7 where n is the dimension of the vector space.</footnote>

<footnote n="8">8 Note that the FEMU method can be prone to noise over-fitting. However, this should not be a concern for our pulmonary model application, as we are not estimating many parameters. Moreover, our method allows us to detect potential noise overfitting. Indeed, for a given noise level, we generate multiple realizations of the random Gaussian noise added to the synthetic measurement and we initialize the optimization process randomly. If after the multiple noise realizations, the estimated values remain reasonably close to the ground-truth, the method does not overfit.</footnote>

<footnote n="14">14 This estimation method is applied in the inverse uncertainty quantification pipeline presented previously. In the scope of this study, the noisy-synthetic-measurements are generated for SNRs of +∞, 20, 10, 5 and 2.5. These choices were based on typical SNR values for CT-scans, ranging from 10 to 20 [Li et al. 2019]. The value +∞ allows us to study the identifiability of parameters in an ideal case, where measurements are not noisy. Additionally, we chose to study the identifiability for a SNR of 5 to model very noisy data, and the value 2.5 corresponds to a “non-physical” case, allowing us to study parametrizations robustness to higher noise levels. The estimated parameters are then compared to the ground-truth. We will call the error of a parameter $\theta_i$ (in percentage) the fraction</footnote>
1

2

3

4

5

6

7

8

9

10

11

Note that classical inverse uncertainty quantification problems involve Bayesian approaches [Nagel et al. 2016; Kerfriden et al. 2019], and sensitivity analysis Sobol approaches [Saltelli et al. 1995; Sobol et al. 2007; Wirthl et al. 2023]. In our case, applying such methods would however be very expensive given our model complexity.

In this work, different parametrizations –i.e. number and type of parameters estimated– are investigated. Our study aims at giving information on the number of parameters that can be identified simultaneously with a good accuracy for future clinical applications, by quantifying the robustness of the chosen parametrizations to measurement and/or noise errors. If the results presented in this article are centered on our specific pulmonary model, our method can be generalized to any model requiring inverse problems.

12

13

14

15

16

17

**Model parameters** Our model includes a substantial number of parameters. As a result, we will very likely not be able to identify all of them. Two sets of parameters are however of particular interest, since they could be used to characterize or predict the evolution of IPF: the pleural pressure, which corresponds in our case to the pleural pressure constant $\tilde{p}_0$ introduced in Equation (33) of Appendix A, and the stiffness, which is represented by the constant $\alpha$ appearing in the skeleton free energy $\Psi_{skel}$ [Laville et al. 2023].

18

19

20

21

22

Throughout this article, we will identify lung stiffnesses, which correspond to constants. However, if $\alpha$ is defined with a constant value, it is not constant throughout the lungs: it is indeed spatially distributed because the energy is scaled by the highly heterogeneous porosity. Note also that in this article, all parameters are identified for synthetic data, in order to validate the proposed estimation pipeline.

23

24

25

26

27

28

29

30

31

32

33

For patients with IPF, different zones appear in the lungs: a healthy zone, which is typically closer to the apex, and fibrosed zones, increasingly stiffer from the apex to the base. In this article, we will investigate different parameterizations, i.e., different number of zones, and hence different values/number of $\alpha$ identified. In practice, to create the different zones, we compute the distance between the apex and the base of the lungs h. For n zones considered, we then define $\Delta y = \frac{h}{n}$, which represents the height of each zone. We then create the n zones of height $\Delta y$. A different value of stiffness $\alpha$ is then associated to each zone. In this article, $\alpha_h$ always refers to the healthy zone. $\alpha_f$ refers to the fibrosed stiffness, when exactly two zones are defined. $\alpha_{f1}$ and $\alpha_{f2}$ then correspond to the highest and the second-highest fibrosed stiffnesses, when at least three zones are defined (i.e., there are at least two fibrosed zones and one healthy zone). $\alpha_{f3}$ corresponds to the smallest fibrosed stiffness, when four zones are defined.

34

35

36

37

38

39

40

41

42

The values of the constants used in the model are presented in Table 1. They were determined from the work of [Patte et al. 2022b], who based her mechanical laws and parameter values on the experiments of [Gibson et al. 1976]. The subscript “e” refers to quantities at end-exhalation, and the subscript “i” refers to quantities at end-inhalation. In the following section, different parametrizations, i.e., different combinations of parameters, are studied in terms of noise and/or model robustness. The chosen parametrization corresponds to the vector $\theta$ previously introduced, i.e., the parameters $\theta$ identified, while the other parameters, fixed at reference values, correspond to $\eta$. Unless otherwise specified, the reference values chosen for $\eta$ correspond to the ground-truth values.

43

## 2.4 Implementation details

44

45

46

47

48

49

**Geometry and mesh** In the scope of this article, we will perform estimation on the lung model developed by [Patte et al. 2022b], with the new boundary conditions introduced in [Peyraut et al. 2024a]. Our model was designed to integrate clinical data, which, as a reminder, consists in one CT-scan at end-exhalation and one CT-scan at end-inhalation. Note, however, that as previously explained, the estimation process is not performed on real clinical data, but on a generic right lung model <fnref n="1" />. This generic model corresponds to a patient at end-exhalation, and was meshed

<footnote n="1">¹https://www.zygote.com</footnote>
<table><caption>Table 1: Parameters used for our computations, based on [Patte et al. 2022b].</caption><thead><tr><th>Name</th><th>Symbol</th><th>Unit</th><th>Value</th></tr></thead><tbody><tr><td rowspan="7">Skeleton energy parameters</td><td>α<sub>h</sub></td><td>[kPa]</td><td>0.09</td></tr><tr><td>α<sub>f</sub></td><td>[kPa]</td><td>0.67</td></tr><tr><td>[ α<sub>h</sub> , α<sub>f</sub> ]</td><td>[kPa]</td><td>[0.09, 0.67]</td></tr><tr><td>[ α<sub>h</sub> , α<sub>f2</sub> , α<sub>f1</sub> ]</td><td>[kPa]</td><td>[0.09, 0.67, 0.9]</td></tr><tr><td>[α<sub>h</sub>, α<sub>f3</sub>, α<sub>f2</sub>, α<sub>f1</sub>]</td><td>[kPa]</td><td>[0.09, 0.5, 0.67, 0.9]</td></tr><tr><td>δ</td><td>[-]</td><td>0.5</td></tr><tr><td>β<sub>1</sub></td><td>[kPa]</td><td>0.6</td></tr><tr><td rowspan="2">Bulk energy parameter</td><td>κ</td><td>[kPa]</td><td>100</td></tr><tr><td>φ<sub>s</sub></td><td>[-]</td><td>0.5</td></tr><tr><td>Pore energy parameter</td><td>ζ</td><td>[-]</td><td>10<sup>-5</sup></td></tr><tr><td>End-exhalation pleural pressure</td><td>P<sub>pl,e</sub></td><td>[kPa]</td><td>-0.5</td></tr><tr><td>End-inhalation pleural pressure</td><td>P<sub>pl,i</sub></td><td>[kPa]</td><td>-2.0</td></tr></tbody></table>

with 32416 tetrahedral elements. The problem solved with our model consists in computing an –unknown– end-inhalation configuration from a –known– end-exhalation configuration.

**Influence of the prone and supine positions on the computation** If we study in a first approach the identifiability of the parameters using two CT-scans as input –since it corresponds to the actual data we have access to–, we can also study the influence on the estimation of the use of additional images. As we would in particular like to investigate the relevance of taking images in prone position as input in addition to images in supine position, the geometry at end-exhalation in prone and supine positions should hence be different. However, we only have access to a generic model at end-exhalation.

In order to create our synthetic measurements, we therefore apply the following steps: (i) we compute the unloaded configuration from the generic end-exhalation configuration, considering zero gravity, i.e., the unloading problem merely corresponds to a scaling, (ii) we compute an end-exhalation configuration in supine position from the unloaded configuration, i.e., we consider that gravity is oriented from the ventral to the dorsal area during the computation, (iii) we compute an end-exhalation configuration in prone position from the unloaded configuration, i.e., we consider that gravity is oriented from the dorsal to the ventral area during the computation, (iv) we compute the end-inhalation configurations in supine and prone positions, from the end-exhalation configurations in supine and prone positions, computed in (ii) and (iii) respectively.

Our input data is composed either of one or two displacement fields. When only data in supine position is considered, there is one displacement field as input –displacement field from end-exhalation to end-inhalation–. When data in both prone and supine positions is considered, two displacement fields –displacement field from end-exhalation to end-inhalation in prone position and in supine position– are used as input.

**Implementation** The problems presented in this article are solved using the finite elements method. The resolution of the finite element problem relies heavily on tools developed in FEniCS [Logg et al. 2012; Alnæs et al. 2015], and a library based on FEniCS, “dolfin_mech” [Genet et al. 2024].

The code developed to produce the figures of this article is currently freely available in interactive Notebooks at: https://apeyraut.gitlabpages.inria.fr/uncertainty-quantification-paper-demos [Peyraut et al. 2024c]. Note that we recommend to launch computations with adequate computation resources. Additionally, the results obtained with the code available online might slightly differ from the results presented in this work. Our results were indeed obtained with the
1 Zygote geometry, which is submitted to copyright and could hence not be freely included in the

2 Notebooks. We proposed a generic cube mesh as a replacement.

3 Results

4 In order to assess the identifiability of our parameters, we assess the robustness of the estimation

5 with respect to noise and to model errors for a simple case, i.e., identifying parameters in a lung

6 with one zone. We then gradually increase the complexity, by identifying more parameters and by

7 increasing the number of zones defined in the lungs.

3.1 Identification errors with one zone

3.1.1 Impact of measurement errors on the identified stiffness

10 We first choose to consider a lung with only one zone, i.e., a healthy patient. We then quantify
11 the identifiability of the stiffness, when all other parameters are fixed at their ground-truth values
12 (detailed in the Methods section). The error distribution of the estimated stiffness as a function
13 of the SNR is presented in Figure 5. In this figure, and for all other figures of this article, the
14 straightlines represent the mean of the distribution, while the transparent areas represent the
15 standard deviations of the distributions.

<figure><img src="image_6.png" /><figcaption>Impact of measurement errors, single-zone lung, identification of $\alpha_h$ only</figcaption></figure>

Figure 5: Error distribution on the estimated stiffness, depending on the SNR, when only one zone—which is healthy— is considered. When no noise is applied to the synthetic data, the stiffness is perfectly identified, and the error committed on the estimated parameter increases as the noise applied to the data increases.

16 Our pipeline allows us to quantify the identifiability of the stiffness depending on the noise present
17 on the measurements for our poromechanical framework, and to quantify the level of confidence
18 we have on our estimated parameter. The results we obtain are classical: the identified stiffness is
19 close to the ground-truth for large SNRs, while the error committed during the estimation process
20 increases for smaller SNRs.

21 In the scope of our study, typical levels of SNRs are around 10 (typical noise values for CT-scans).
22 From our results, it is therefore possible to access to a confidence interval and to anticipate that
23 for a lung with one zone, the stiffness can be estimated with a confidence of about 2 %. This
24 information is of particular interest for clinical applications.

3.1.2 Impact of measurement errors on the identified pleural pressure

26 Since the pleural pressure is also a very interesting parameter to identify in our problem, we also

27 quantify its noise robustness. To do so, all parameters, except the end-inhalation pleural pressure

28 $p_{pl,i}$, are fixed at reference values, and $p_{pl,i}$ is then estimated for different noise levels. The obtained

29 results are presented in Figure 6.
Impact of measurement errors, single-zone lung, identification of $p_{pl,i}$ only

<figure><img src="image_7.png" /><figcaption>Figure 6: Error distribution on the estimated end-inhalation pleural pressure $p_{pl,i}$, for a lung with one zone, depending on the SNR. The error committed during the estimation is very low, for all SNRs investigated, even though it increases slightly for smaller SNRs.</figcaption></figure>

The error committed on the estimated end-inhalation pleural pressure is very low, even for small SNRs. For the highest noise level -SNR=2.5-, the parameter is indeed identified with a confidence of about +/- 0.5%. The pressure identification seems to be slightly more robust than the stiffness identification, as the error seems to be less noise sensitive. Additionally, both the stiffness and the end-inhalation pleural pressure seem to be identifiable, when all other parameters are fixed at their ground-truth values. As a result, we can now try to identify both parameters simultaneously, while fixing all other parameters are reference values.

### 3.1.3 Simultaneous estimation of the stiffness and of the pleural pressure

The simultaneous estimation of $\alpha$ and $p_{pl,i}$ for a lung with one zone is presented in the top row of Figure 7. As there is only one zone, we choose the stiffness as healthy ($\alpha_0 = \alpha_h$). Identifying both parameters simultaneously seems to give a good estimation quality for large SNRs.

Both parameters are identifiable: when no noise is applied to the synthetic measurements, both the stiffness $\alpha_h$ and the pleural pressure $p_{pl,i}$ are perfectly identified. For a relatively large SNR, e.g., 20, the error remains small: under 6% for $\alpha_h$ and under 2% for $p_{pl,i}$. However, the error committed for typical noise levels, i.e., a SNR of 10, is much larger: an error of about 15% is observed for $\alpha_h$ and 5% for $p_{pl,i}$.

These high errors for a typical CT-scan SNR value might be explained by the low identifiability of our problem, as we are trying to identify both the stiffness and the loading. In general, both parameters should not be simultaneously identifiable. However, the high heterogeneity of the lungs, due to the highly heterogeneous porosity, leads to a highly inhomogeneous stiffness in the lungs, which explains why we can identify simultaneously both parameters. We can yet suppose that this identification is not very robust, which is what is indeed observed when noise increases.

Note that it is impossible to identify, for example, both pleural pressures -end-exhalation and end-inhalation- simultaneously. Only the differential in pleural pressure would indeed be identified. Additionally, our model also seems to display different sensitivities to the different parameters. In particular, given the results obtained previously, we could expect a higher sensitivity of the model to the pleural pressure than to the stiffness.

As the model seems sensitive to $p_{pl,i}$, we can therefore expect a high impact on the estimation of model errors on $p_{pl,e}$. To investigate this hypothesis, let us therefore study the impact of model errors on $p_{pl,e}$ on estimation errors on $\alpha_h$ and $p_{pl,i}$. To do so, the estimation process is performed when the end-exhalation pleural pressure $p_{pl,e}$ is fixed at -20, -10, 0, 10, 20 % of its ground-truth value during the estimation process.
<figure><img src="image_8.png" /><figcaption>Impact of model errors, single-zone lung, identification of ($\alpha_h, p_{pl,i}$)</figcaption></figure>

<figure><img src="image_9.png" /></figure>

Figure 7: Error distribution for the estimated stiffness (left) and the end-inhalation pleural pressure (right), depending on the SNR (top row), when both parameters are identified at the same time. When no noise is applied to the synthetic data, both parameters are perfectly identified. The error committed on the estimated parameters increases as the noise applied to the data increases. The error on the stiffness seems to be higher than the error on the end-inhalation pleural pressure. When a bias is introduced on the fixed end-exhalation pleural pressure (bottom row) the mean estimation error for a SNR of 10 increases significantly for both parameters with errors on $p_{pl,e}$.

The results are generated for a SNR of 10 –typical noise levels for CT-scans– and are presented in the bottom row of Figure 7. As expected, our model displays a high sensitivity to errors committed on $p_{pl,e}$. For example, for an error of 10 % on $p_{pl,e}$, the average error is about 4% for $\alpha$, and about 15% for $p_{pl,i}$. The error committed on $p_{pl,i}$ is globally higher than for $\alpha$. Additionally, an overestimation of $p_{pl,e}$ leads to an increase of the error approximately equivalent –even though slightly higher– for the estimated $p_{pl,i}$. For example, overestimating $p_{pl,e}$ by 20% leads to an overestimation of $p_{pl,i}$ of about 30%.

As a result, this parametrization displays a high noise sensitivity, as the error committed during the estimation process is quite high for SNRs of 10. Additionally, our model seems to be quite sensitive to model errors on $p_{pl,e}$. As it is impossible to identify both pleural pressures, and as the end-exhalation pleural pressure might vary from a patient to another, this represents a limitation of our model that we should be aware of during the estimation process for clinical applications.
1

## 3.2 Identification errors with two zones

2

3

4

5

6

7

If identifying both the loading and the stiffness simultaneously does not give good quality results, the estimation quality is however very good when the parameters are identified separately. As a result, we can try to perform the estimation when the lung is separated into two zones –a fibrosed zone of stiffness $\alpha_f$ at the base of the lungs, and a healthy zone at the apex, with a stiffness $\alpha_h$, if both stiffnesses of $p_{pl,i}$ are not estimated simultaneously. Fixing one of the three parameters at its ground-truth value should indeed increase the identifiability of the problem.

8

9

10

11

12

13

14

Different parametrizations have been investigated in the past. Notably, in [Patte et al. 2022a], both stiffnesses are identified, while both pleural pressures are fixed at a physiological values. In [Laville et al. 2023], the fibrosed stiffness and the end-inhalation pleural pressure are identified, while all other parameters are fixed at physiological values. When the estimation was performed, we, however, did not have a clear idea on the actual identifiability and the reliability of the estimated parameters. In order to validate our previous parametrizations, we therefore decided to apply our inverse uncertainty quantification pipeline to both cases.

15

### 3.2.1 Analysis of [Patte et al. 2022a]'s parametrization

16

17

18

19

20

For this parametrization, both pleural pressures $-p_{pl,i}$ and $p_{pl,e}$ are fixed at reference values, chosen within a physiological range. The healthy stiffness $\alpha_h$ and the fibrosed stiffness $\alpha_f$ are then simultaneously identified. In order to quantify the reliability of this parametrization, we identify simultaneously the healthy and fibrosed stiffnesses for different noise levels. The results are presented in the first row of Figure 8.

21

22

23

24

When all non-identified parameters are fixed at their ground-truth value, the estimation quality is good for both identified stiffnesses. For a SNR of 10, the healthy stiffness is identified perfectly, with an accuracy of +/−4%, while the fibrosed is identified perfectly as well, with an accuracy of +/−3%. The standard deviation of the error distributions increases with noise.

25

26

27

28

29

30

31

We also investigate the impact of model errors on the estimation. In particular, as both pleural pressures are fixed at reference values, we investigate the impact of errors on $p_{pl,i}$ on the estimation quality. To do so, both stiffnesses are simultaneously identified for a SNR of 10, when an error of -20 %, -10 %, 0 %, 10 %, and 20 % is committed on $p_{pl,i}$ during the estimation process. The obtained results are presented in the bottom row of Figure 8. Errors on the end-inhalation pleural pressure have a significant impact on the estimation of both parameters. For example, an error of +10% on $p_{pl,i}$ leads to an error of 37% +/−5% for $\alpha_h$ and of 28% +/−9% for $\alpha_f$.

32

33

34

35

36

37

38

Both stiffnesses behave similarly: when the $p_{pl,i}$ is overestimated (i.e., it is more negative), both stiffnesses increase. This result seems consistent: if the pleural pressure becomes more negative, and if all other parameters are fixed, the displacement from end-exhalation to end-inhalation becomes higher. As a result, to retrieve the synthetic measured displacement field, the stiffnesses of the lung should increase. If both stiffnesses display the same overall behavior, $\alpha_h$ seems to be a little more sensitive to model errors, as the average of the distribution is consistently further from the ground-truth than for $\alpha_f$.

39

40

41

If the parametrization proposed by [Patte et al. 2022a] leads to a good noise robustness, it is however very sensitive to model errors committed on the fixed pressure. As a result, this parametrization may not be optimal.

42

### 3.2.2 Analysis of [Laville et al. 2023]'s parametrization

43

44

45

46

47

48

49

For the proposed parametrization, the fibrosed stiffness $\alpha_f$ and the end-inhalation pleural pressure $p_{pl,i}$ are identified simultaneously, while fixing all other parameters, including the healthy stiffness $\alpha_h$ at physiological values. The error distributions for the identified parameters $\alpha_f$ and $p_{pl,i}$ depending on the noise are presented in Figure 9, in the first row. The error committed on the identified $p_{pl,i}$ is very small, while it is higher for the estimated $\alpha_f$ for higher noise levels. However, for all noise levels, the mean of the distribution is very close to the ground-truth for both estimated parameters. Additionally, the error for the identified $\alpha_f$ remains fairly small for all SNRs. For
Impact of measurement errors, two-zone lung,
parametrization from [Patte et al. 2022a]: identification of ($\alpha_h, \alpha_f$)

<figure><img src="image_10.png" /></figure>

Impact of model errors, two-zone lung,
parametrization from [Patte et al. 2022a]: identification of ($\alpha_h, \alpha_f$)

<figure><img src="image_11.png" /><figcaption>Figure 8: Error distribution for the simultaneously estimated healthy (left) and fibrosed (right) stiffness, depending on the noise level (top row) and depending on the error on $p_{pl,i}$ (bottom row) for a SNR of 10. For standard CT-scans levels of noise, both parameters can be identified with a good accuracy. They are however very sensitive to model errors.</figcaption></figure>

1. example, for a SNR of 2.5 the fibrosed stiffness is identified perfectly, with an uncertainty of about 15%.

2. For this study, $p_{pl,e}$ was fixed at its ground-truth value. We, however, previously showed that our model seems sensitive to model errors on the pleural pressure. Additionally, as one of the two stiffnesses is fixed, we can investigate the impact of errors on $\alpha_h$ and on $p_{pl,e}$ on the estimation of $\alpha_f$ and of $p_{pl,i}$. To do so, we fix $\alpha_h$, at -20 %, -10 %, 0 %, +10 % and +20 % of its ground truth value during the estimation, while fixing all other parameters at their reference values. The same study is also conducted, by fixing the end-exhalation pleural pressure $p_{pl,e}$ at a wrong value. The obtained results are presented in the bottom rows of Figure 9.

3. $p_{pl,i}$ displays a very good robustness to model errors on the fixed healthy stiffness. For example, for an error of +10% on the healthy stiffness, the error committed on the identified $p_{pl,i}$ is around 1.6% +/- 0.8%. $\alpha_f$ seems, however, much more sensitive to errors on $\alpha_h$: an error of +10% on the healthy stiffness indeed leads to an error of about 7% +/- 4% on the estimated $\alpha_f$.

4. Our estimation additionally seems to remain very sensitive to errors on $p_{pl,e}$. For an error of +10% on $p_{pl,e}$, $\alpha_f$ is estimated with a relatively small error (around 0% +/- 5%). The error committed
Impact of measurement errors, two-zone lung,
parametrization from [Laville et al. 2023]: identification of ($\alpha_f$, $p_{pl,i}$)

<figure><img src="image_12.png" /></figure>

Impact of model errors, two-zone lung,
parametrization from [Laville et al. 2023]: identification of ($\alpha_f$, $p_{pl,i}$)

<figure><img src="image_13.png" /><figcaption>Figure 9: Error distribution for the simultaneously estimated fibrosed stiffness (left) and end-inhalation pleural pressure (right), depending on the noise level (top row), and on model errors on $\alpha_h$ (middle row) and on $p_{pl,e}$ (bottom row), for a SNR of 10. For standard CT-scans levels of noise, both parameters can be identified with a good accuracy when all other parameters are fixed at their ground-truth value. Additionally, our model displays a higher sensitivity to errors on $p_{pl,e}$ than to errors on $\alpha_h$.</figcaption></figure>

<footnote n="1">1 on $p_{pl,i}$ is however much more significant: around 16% +/- 4%. Similarly to what was observed for a lung with one zone, $p_{pl,i}$ is systematically slightly more misestimated than $p_{pl,e}$.</footnote>

<footnote n="2">2 As the error committed on the $\alpha$ is already high when both the $\alpha$ and $p_{pl,i}$ are simultaneously identified for a one-zone lung, it is not relevant to try to identify all three parameters $\alpha_h$, $\alpha_f$ and $p_{pl,i}$ simultaneously. We can, however, attempt to identify more parameters, by increasing the number of zones and fixing $\alpha_h$ at a reference value to increase the identifiability of the problem.</footnote>
## 3.3 Estimation with n zones

### 3.3.1 Estimation for a lung with 3 zones

Let us define three zones in the lung. The zone closer to the base of the lung is associated with the highest stiffness, $\alpha_{f1}$. The fibrosed stiffness $\alpha_{f2}$ is associated with the mid-zone, while the zone closer to the apex corresponds to the healthy zone $\alpha_h$. As mentioned earlier, the estimation of all stiffnesses and of $p_{pl,i}$ will not lead to a good estimation quality. As a result, we can attempt first to identify all fibrosed stiffnesses $-\alpha_{f1}$ and $\alpha_{f2}$ along with $p_{pl,i}$, with all other parameters, including $\alpha_h$, fixed at reference values. The error distributions for each parameter depending on the noise are presented in Figure 10.

<figure><img src="image_14.png" /><figcaption>Figure 10: Error distribution for the simultaneously estimated parameters $\alpha_{f1}$ (top left), fibrosed stiffness $\alpha_{f2}$ (top right), and $p_{pl,i}$ (bottom), depending on the noise level. The pleural pressure is still identified with good accuracy. However, the stiffnesses are not identifiable anymore.</figcaption></figure>

The identification of all fibrosed stiffnesses and of the end-inhalation pleural pressure simultaneously does not give good results. In particular, an initial uncertainty appears for the stiffnesses when no noise is added to the data: +/-4% for $\alpha_{f1}$, +/-8% for $\alpha_{f2}$, and +/-1% for $p_{pl,i}$. For typical noise levels for CT-scans, the estimation remains, however, robust. For a SNR of 10, the stiffnesses are identified roughly with the same uncertainty as for a SNR of +∞ even though our results suggest that the stiffnesses are not uniquely identifiable. As the error for a SNR of 10 can still be deemed acceptable for all parameters (as it remains under 10%), we can try to perform the estimation for a lung with more zones.
### 3.3.2 Estimation for a lung with 4 zones

In this case, the lung is separated in 4 zones: a healthy zone at the apex, and three fibrosed zones of stiffness $\alpha_{f1}$ at the base of the lungs, $\alpha_{f2}$ in adjacent zone and $\alpha_{f3}$ in the zone adjacent to the healthy zone, such that $\alpha_{f3} < \alpha_{f2} < \alpha_{f1}$. We then identify all fibrosed stiffnesses, along with $p_{pl,i}$, while fixing all other parameters, including $\alpha_h$, at reference values. The obtained results are presented in Figure 11.

<figure><img src="image_15.png" /><figcaption>Impact of measurement errors, four zone-lung, estimation of ($\alpha_{f1}$, $\alpha_{f2}$, $\alpha_{f3}$, $p_{pl,i}$)</figcaption></figure>

Figure 11: Error distribution for the simultaneously estimated fibrosed stiffness $\alpha_{f1}$ (top left), fibrosed stiffness $\alpha_{f2}$ (top right), fibrosed stiffness $\alpha_{f3}$ (bottom left) and end-inhalation pleural pressure (bottom right), depending on the noise for images in supine position for a lung with four zones. The quality of the estimation is poor. In particular, an uncertainty of almost +/ - 30% appears for $\alpha_{f3}$ for all noise levels.

The estimation quality is not good for all parameters except for $p_{pl,i}$. The standard deviation of the error distribution is approximately the same as for the initialization for $\alpha_{f3}$ for all noise levels; the parameter is not identifiable. The standard deviation of the distributions is smaller for the other stiffnesses, and remains roughly the same for all SNRs: about 12% for $\alpha_{f1}$, and 15% for $\alpha_{f2}$. The error committed on $p_{pl,i}$, however, remains very small. As the stiffnesses are not identifiable, since they are not uniquely identified in the absence of noise, and as the standard deviation is very high for $\alpha_{f3}$, it is not relevant to increase the number of zones studied, as it will lead to more uncertainty on the identified parameters.

## 3.4 Improving the estimation process using gravity

As a result, the parametrization proposed by [Laville et al. 2023] seems the most relevant when high accuracy is required and when the only data available consists in one CT-scan at end-exhalation
and one at end-inhalation for a patient in supine position. However, we could potentially access to images in prone position as well. Incorporating more clinical data during the estimation, e.g., the displacement field from end-exhalation to end-inhalation in prone position and in supine positions, could hopefully improve this process.

In order to assess whether including images in prone position in our dataset would improve the estimation, we modified slightly our cost function. If it initially solely included a displacement field from end-exhalation to end-inhalation in supine position, the modified cost function takes additionally in input the displacement field from end-exhalation to end-inhalation in prone position.

### 3.4.1 Estimation of all stiffnesses and of the pleural pressure

We previously managed to identify both $p_{pl,i}$ and $\alpha$ simultaneously, when only one zone was considered. However, this parametrization did not display a very good noise robustness. As a result, we did not attempt to perform the estimation on a more complex case: a lung with two zones –with a fibrosed stiffness $\alpha_f$ at the base and a healthy stiffness $\alpha_h$ at the apex–, where both $\alpha_h$ and $\alpha_f$ would be identified along with $p_{pl,i}$.

However, since the proposed method to use more images to improve the estimation quality may increase the identifiability of our parameters, we apply it on this case. As a reminder, even if it is counterintuitive to expect the identifiability of the three parameters simultaneously, the high heterogeneity of the material may allow us to identify them all. The error distributions for the identified parameters depending on the noise are presented in Figure 12. The distribution in blue, in red and in green are respectively associated with results obtained for a patient in supine position only, in prone position only, and both in prone and supine positions.

For a patient in prone or supine position only, the three parameters are identifiable, since the standard deviation of the error distribution, as well as the average of the distribution are close to 0 when no noise is added to the synthetic data. However, as noise increases, the error committed on the estimated parameters becomes high. For example, for a SNR of 10, for a patient in supine position, $\alpha_h$ is identified with an error of 24% +/- 80%, $\alpha_f$ with an error of 22% +/- 60% and $p_{pl,i}$ with an error of 0.8% +/- 15%. As a result, the estimation is not robust.

However, combining images in prone and supine position improves drastically the estimation quality. With this dataset, for a SNR of 10, $\alpha_h$ is identified with an error of 2% +/- 7%, $\alpha_f$ with an error of 1% +/- 2% and $p_{pl,i}$ with an error of 0.5% +/- 1%. These errors might be acceptable for clinical applications. Enriching our dataset therefore leads to a drastic improvement of the estimation process, by making more parameters identifiable.

### 3.4.2 Estimation with n zones

Although this new method improves the estimation, it does not allow us to increase the number of zones for which the parameters are identifiable. For example, Figure 13 presents the error distribution depending on the noise of the identified parameters $\alpha_{f1}$, $\alpha_{f2}$ and $p_{pl,i}$ for a three-zone lung when all other parameters are fixed at reference values. The obtained results are very close to the case where the only data available is from a patient in supine position. Using images in both prone and supine position does therefore not improve the estimation quality in this case. As the results are very close to the case using images in supine position only, it is not relevant to study the impact on the estimation quality of enriching the dataset for a lung with more zones.

As a result, enriching the dataset allows us to improve the identification. In particular, parametrizations that are not identifiable with images in supine position only may become identifiable when images both in prone and in supine positions are included in the cost function. It is the case, for example, of a lung with two zones, when both stiffnesses and the end-inhalation pleural pressure are simultaneously identified. Enriching the dataset with images in prone position does, however, not allow us to increase greatly the number of identifiable parameters. Our results indeed suggest that for more than three zones, identifying all fibrosed stiffnesses and the end-inhalation pleural pressure will lead to a very poor estimation quality. The improvements yield by this method –e.g.,
Impact of measurement errors, two-zone lung,
new parametrization: identification of ($\alpha_h, \alpha_f, p_{pl,i}$)

<figure><img src="image_16.png" /><figcaption>Figure 12: Error distribution for the simultaneously estimated healthy stiffness $\alpha_h$ (top left), fibrosed stiffness $\alpha_f$ (top right) and end-inhalation pleural pressure $p_{pl,i}$ (bottom), depending on the noise level. All three parameters are identifiable, for all input datasets considered. The error is however very high for a SNR of 10 for a patient in prone and supine positions only. Combining information from the two positions however seems to lead to a very good parameter identifiability for all noise levels.</figcaption></figure>

1. identifying an additional parameter for a lung with two zones— remain however interesting, and
2. allow us to improve the overall estimation quality.

## 4 Discussion

The proposed pipeline investigates the identifiability of the parameters of our model by quantifying the robustness to noise and model errors of different parametrizations, enabling the identification of the optimal set of parameters to identify, and providing confidence intervals for the estimated parameters. This method allows us to study the actual parameter identifiability when using noisy data, which is particularly interesting for biomechanical applications [Marquis et al. 2018]. Applied to a poromechanical lung model, the method improves reliability for clinical applications, by providing insights on what parameters to identify. In our case, the best parametrization corresponds to the simultaneous identification of $p_{pl,i}$ and $\alpha_f$ for a two-zone lung.

12
<figure><img src="image_17.png" /><figcaption>Figure 13: Error distribution for the simultaneously estimated fibrosed stiffness $\alpha_{f1}$ (top left), fibrosed stiffness $\alpha_{f2}$ (top right) and end-inhalation pleural pressure (bottom), depending on the noise for images in supine position (blue) and images in prone and supine position (green) for a lung with three-zone. Enriching the dataset does not seem to improve the estimation process.</figcaption></figure>

# 4.1 Identifiability and sensitivity analysis

## 4.1.1 Parameter identifiability

In the scope of our study, two parameters are of particular interest, as biomarkers could be derived from them: the end-inhalation pleural pressure $p_{pl,i}$ and the stiffness $\alpha$ of the different zones of the lung. We first confirmed that in the most simple case, i.e., a lung with one zone, both the end-inhalation pleural pressure and the stiffness could be identified separately with a good accuracy for typical noise levels, when all other parameters are fixed at their ground-truth value. However, when both parameters are simultaneously identified, the robustness to noise deteriorates: if both parameters are identifiable for reasonable noise levels, the error committed on $\alpha_h$ and $p_{pl,i}$ increases rapidly with noise.

The observed low noise robustness is not surprising in this case. Being able to identify simultaneously both the stiffness $\alpha$ and the end-inhalation pleural pressure $p_{pl,i}$ may, however, seem counterintuitive at first. For a homogeneous material, e.g., for classical hyperelastic problems, it would indeed be impossible to identify both the stiffness and the loading: if the loading is overestimated, underestimating the stiffness should allow us to retrieve a displacement field close to the measurement, and conversely. The pair ($\alpha$, $p_{pl,i}$) is therefore non-unique, and the parameters should not be simultaneously identifiable.
In our poromechanical problem, although stiffness is presented as a constant, it is in practice highly heterogeneous, as it is multiplied by the –highly heterogeneous– porosity, which evolves with applied loading. Misestimating the loading will result in a porosity field different from the ground-truth. Contrarily to heterogeneous material, simply adjusting $\alpha$ will very unlikely lead to a displacement field close to the measurement, due to the different porosity field. As a result, there should be a unique ($\alpha, p_{pl,i}$) matching the –synthetic– measured displacement field, which explains the simultaneous identifiability of both parameters. However, since the identifiability relies solely on the porosity, the estimation is not be very robust, and the error increases rapidly with noise.

As the estimation quality was good for a lung with one zone, we investigated the identifiability of more complex parametrizations, by increasing the number of zones in the lung. However, based on previous observations, all stiffnesses and the loading will not be simultaneously identified, due to low noise robustness. For all parametrizations, we therefore choose to fix $\alpha_h$ at a reference value, and identify $p_{pl,i}$ and the fibrosed stiffnesses, as we are more interested in the evolution of diseased tissues.

The obtained results show that our problem is identifiable for a lung with up to two zones. For a lung with three zones, the problem is indeed ill-posed, as there is not a unique solution for stiffnesses when no noise is added to the synthetic measure. The error on the estimated parameters remains, however, stable for all noise levels, suggesting that the local minima are very close to the ground-truth value. Depending on the accuracy level required for the estimated parameters, this parametrization may remain an acceptable option, as the error on the estimated parameters remains under 15% for typical noise levels.

However, for a lung with four zones, it is impossible to identify all fibrosed stiffnesses and $p_{pl,i}$ simultaneously. If the pleural pressure is still identified with a good accuracy, and if the error obtained for the two stiffer stiffnesses might be considered acceptable, the error on the smaller stiffness is very high: it is close to the initialization interval for all SNRs. As a result, this parametrization does not give good quality results.

As mentioned in the Methods section, the FEMU method, used for the estimation in this work, is prone to overfitting. However, based on our results, we can conclude that the method does not overfit when applied to our pulmonary model, when the lung has one or two zones but tends to overfit when more zones are considered. For these cases, additional a priori information should be included. While regularization techniques could be considered to address this issue [Ahmadian et al. 1998; Friswell et al. 2001; Hua et al. 2009], exploring them is beyond the scope of this article.

The identification problem for a lung with three zones or more is ill-posed, as the parameters cannot be uniquely identified with zero noise on the data [Ghattas et al. 2021]. As a result, a Bayesian approach could provide interesting results for these cases, as they would provide a probabilistic distribution on the estimated parameters and would make the estimation more robust [Santos et al. 2009]. However, these approaches require the knowledge of a priori parameters distributions [Wu et al. 2019], which are unknown in our case, and could bias the estimation.

### 4.1.2 Sensitivity analysis

Our model displays different sensitivities to the parameters. In particular, our model seems more sensitive to the pleural pressure than to the stiffness. Errors on the fixed $p_{pl,e}$ indeed lead to higher estimation errors on the estimated parameters than errors on $\alpha$, as can be seen in Figure 9, for example. Additionally, errors committed on the estimated end-inhalation pleural pressure due to noise are globally smaller than on the estimated stiffness. This observation is aligned with [Wu et al. 2019], who studied the link between sensitivity and identifiability, and concluded that there seems to be a close link between identifiability and high sensitivity.

Moreover, the stiffnesses seem to be identified with a higher accuracy for higher stiffnesses than for smaller stiffnesses, which is noticeable in particular for the case of the lungs with four zones, for which the error on the two highest stiffnesses remains below 15%, while it is around 30% for the lowest stiffness.
1

2

2015], who concludes that uncertainties on their model parameters have a greater impact on the estimation than uncertainties on measurements.

3

## 4.2 Improving the estimation process by enriching the input dataset

4

5

6

7

8

9

10

The results previously discussed are obtained by using data from a generic patient in supine position only, which corresponds the clinical data we have access to for our clinical studies. However, enriching the dataset by also using data from a patient in prone position allows us to improve the estimation quality. As the displacement field differs for both positions, using images in both positions indeed allows us to suppress potential minima: a minimum in supine position is very unlikely a minimum in prone position, and conversely. Our method therefore allows us to increase the identifiability of our parameters.

11

12

13

This is what is indeed observed for a lung with two zones, when $\alpha_h$, $\alpha_f$ and $p_{pl,i}$ are simultaneously identified. When using images in supine position only, the problem is not identifiable. However, when combining images in both positions, the problem becomes identifiable.

14

15

16

17

18

It is however important to note that this method is only efficient because the impact of gravity is noticeable on the lung behavior. If the loading was too high compared to the effect of gravity, using images in prone and in supine positions would not have an effect on the estimation quality, as there would not be noticeable differences between the two cases. This method seems however very promising for improving the estimation quality in the case of the lungs.

19

# 5 Conclusion

20

21

22

23

24

25

26

Parameter estimation is crucial for many model personalization problems, but its quality depends strongly on the available data. Furthermore, it is impossible to quantify the error in estimated parameters without prior studies. To address this issue, we propose in this article a pipeline to quantify the robustness to noise and model errors for specific parameterizations. Our method identifies the best parametrization for a problem, quantifies model sensitivity to parameters, and provides a confidence interval for the estimated parameters. Our pipeline was applied to a poromechanical lung model, but could be extended to any model.

27

28

29

30

31

32

33

Our method relies on the creation of synthetic noisy measurements and performs the estimation, potentially biased by model errors, using the FEMU method. Our results provide information on the reliability of different parametrizations, which is particularly valuable for clinical applications. For our *in-silico* lung model, the reliability of the estimation could not be determined without a prior study. Our study provides access to quantitative information allowing us to choose the most relevant parametrization, based on the confidence intervals identified for different noise levels and the robustness to model errors.

34

35

36

37

38

39

40

All results presented are obtained with the FEMU method, which is computationally expensive, due to the model evaluation at each step of the optimization. Alternative methods, such as direct methods like the Equilibrium Gap Method (EGM), or methods at the interface of identification and tracking [Peyraut et al. 2025], could be considered. They are however not implemented in our problem yet, since our starting position is already loaded; applying direct methods is therefore not straightforward. Model reduction could also be considered to reduce computation costs [Barahona et al. 2024].

41

42

43

44

45

46

47

48

Our study could also be extended using more formal methods, such as Bayesian approaches or Sobol analyses. However, if these approaches apply well to simple models, e.g., problems in linear elasticity, they are expensive for non-linear models. We could work on a surrogate model [Tosin et al. 2020] in order to apply these techniques. Bayesian approaches using residuals for our equations could also be studied to make the estimation less expensive [Cheung et al. 2019]. Additionally, a study on applying the EGM to prestressed structure could help speeding the identification process. In the meantime, the proposed method provides key information on parameter estimation for complex problems.

49

Our pipeline can therefore be used to improve the estimation process, by studying the identifiability
4 Acknowledgements

<footnote n="1">1 of the different parameters for different modeling choices, model parametrizations, identification</footnote>

<footnote n="2">2 approaches, etc. This pipeline is therefore an important step towards improving the reliability of</footnote>

<footnote n="3">3 our personalized lung model.</footnote>

<footnote n="5">5 We would like to thank the French National Research Agency (ANR) for providing financial support</footnote>

<footnote n="6">6 to this project (Grant no: ANR-19-CE45-0007).</footnote>
References
1

2 Aguado-Sierra, J., A. Krishnamurthy, C. Villongco, J. Chuang, E. Howard, M. J. Gonzales, J. Omens, D. E. Krummen, S. Narayan, R. C. Kerckhoffs, and A. D. McCulloch (2011). "Patient-Specific Modeling of Dyssynchronous Heart Failure: A Case Study". In: *Progress in Biophysics and Molecular Biology* 107.1, pp. 147-155. DOI: 10.1016/j.pbiomolbio.2011.06.014.

3 Agustí, A., E. Melén, D. L. DeMeo, R. Breyer-Kohansal, and R. Faner (2022). "Pathogenesis of Chronic Obstructive Pulmonary Disease: Understanding the Contributions of Gene-Environment Interactions across the Lifespan". In: *The Lancet Respiratory Medicine* 10.5, pp. 512-524. DOI: 10.1016/S2213-2600(21)00555-5.

10 Ahmadian, H., J. Mottershead, and M. Friswell (1998). "Regularisation Methods For Finte Element Model Updating". In: *Mechanical Systems and Signal Processing* 12.1, pp. 47-64. DOI: 10.1006/mssp.1996.0133.

13 Alnæs, M., J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, and G. N. Wells (2015). "The FEniCS Project Version 1.5". In: *<p>Archive of Numerical Software* Vol 3. DOI: 10.11588/ANS.2015.100.20553.

16 Anstett-Collin, F., L. Denis-Vidal, and G. Millérioux (2020). "A Priori Identifiability: An Overview on Definitions and Approaches". In: *Annual Reviews in Control* 50, pp. 139-149. DOI: 10.1016/j.arcontrol.2020.10.006.

19 Antoniou, K. M., G. A. Margaritopoulos, S. Tomassetti, F. Bonella, U. Costabel, and V. Poletti (2014). "Interstitial Lung Disease". In: *European Respiratory Review* 23.131, pp. 40-54. DOI: 10.1183/09059180.00009113.

22 Avril, S., M. Grédiac, and F. Pierron (2004). "Sensitivity of the Virtual Fields Method to Noisy Data". In: *Computational Mechanics* 34.6, pp. 439-452. DOI: 10.1007/s00466-004-0589-6.

24 Avril, S., M. Bonnet, A.-S. Bretelle, M. Grédiac, F. Hild, P. Ienny, F. Latourte, D. Lemosse, S. Pagano, E. Pagnacco, and F. Pierron (2008). "Overview of Identification Methods of Mechanical Parameters Based on Full-field Measurements". In: *Experimental Mechanics* 48.4, pp. 381-402. DOI: 10.1007/s11340-008-9148-y.

28 Barahona, J., F. S. Costabal, and D. E. Hurtado (2024). "Machine learning modeling of lung mechanics: Assessing the variability and propagation of uncertainty in respiratory-system compliance and airway resistance". In: *Computer Methods and Programs in Biomedicine* 243, p. 107888. DOI: 10.1016/j.cmpb.2023.107888.

32 Barrau, N. (2024). "3D MR Spirometry". PhD thesis. Université Paris-Saclay.

33 Bates, J. H. T., G. S. Davis, A. Majumdar, K. J. Butnor, and B. Suki (2007). "Linking Parenchymal Disease Progression to Changes in Lung Mechanical Function by Percolation". In: *American Journal of Respiratory and Critical Care Medicine* 176.6, pp. 617-623. DOI: 10.1164/rccm.200611-17390C.

37 Bordas, R., C. Lefevre, B. Veeckmans, J. Pitt-Francis, C. Fetita, C. E. Brightling, D. Kay, S. Siddiqui, and K. S. Burrowes (2015). "Development and Analysis of Patient-Based Complete Conducting Airways Models". In: *PLOS ONE* 10.12. Ed. by J. Sznitman, e0144105. DOI: 10.1371/journal.pone.0144105.

41 Boucneau, T., B. Fernandez, P. Larson, L. Darrasse, and X. Maître (2020). "3D Magnetic Resonance Spirometry". In: *Scientific Reports* 10.1, p. 9649. DOI: 10.1038/s41598-020-66202-7.

43 Carniel, E. L., I. Toniolo, and C. G. Fontanella (2020). "Computational Biomechanics: In-Silico Tools for the Investigation of Surgical Procedures and Devices". In: *Bioengineering* 7.2, p. 48.

45 Caruel, M., R. Chabiniok, P. Moireau, Y. Lecarpentier, and D. Chapelle (2014). "Dimensional Reductions of a Cardiac Model for Effective Validation and Calibration". In: *Biomechanics and modeling in mechanobiology* 13, pp. 897-914.

48 Chabiniok, R., P. Moireau, P.-F. Lesault, A. Rahmouni, J.-F. Deux, and D. Chapelle (2012). "Estimation of Tissue Contractility from Cardiac Cine-MRI Using a Biomechanical Heart Model". In: *Biomechanics and modeling in mechanobiology* 11, pp. 609-630.

51 Chapelle, D. and P. Moireau (2014). "General Coupling of Porous Flows and Hyperelastic Formulations—from Thermodynamics Principles to Energy Balance and Compatible Time Schemes". In: *European Journal of Mechanics-B/Fluids* 46, pp. 82-96.

54 Cheung, S. W. and N. Guha (2019). "Dynamic Data-Driven Bayesian GMsFEM". In: *Journal of Computational and Applied Mathematics* 353, pp. 72-85. DOI: 10.1016/j.cam.2018.12.010.
1

Claire, D., F. Hild, and S. Roux (2004). "A Finite Element Formulation to Identify Damage Fields: The Equilibrium Gap Method". In: *International Journal for Numerical Methods in Engineering* 61.2, pp. 189-208. DOI: 10.1002/nme.1057.

2

3

Cook, D., M. Julias, and E. Nauman (2014). "Biological Variability in Biomechanical Engineering Research: Significance and Meta-Analysis of Current Modeling Practices". In: *Journal of Biomechanics* 47.6, pp. 1241-1250. DOI: 10.1016/j.jbiomech.2014.01.040.

4

5

6

Coorey, G., G. A. Figtree, D. F. Fletcher, V. J. Snelson, S. T. Vernon, D. Winlaw, S. M. Grieve, A. McEwan, J. Y. H. Yang, P. Qian, K. O'Brien, J. Orchard, J. Kim, S. Patel, and J. Redfern (2022). "The Health Digital Twin to Tackle Cardiovascular Disease—a Review of an Emerging Interdisciplinary Field". In: *npj Digital Medicine* 5.1, p. 126. DOI: 10.1038/s41746-022-00640-7.

7

8

9

10

Coussy, O. (2004). *Poromechanics*. John Wiley & Sons.

12

13

Daly, A. C., D. Gavaghan, J. Cooper, and S. Tavener (2018). "Inference-Based Assessment of Parameter Identifiability in Nonlinear Biological Models". In: *Journal of The Royal Society Interface* 15.144, p. 20180318. DOI: 10.1098/rsif.2018.0318.

14

15

Demiray, H. (1972). "A Note on the Elasticity of Soft Biological Tissues". In: *Journal of biomechanics* 5.3, pp. 309-311.

16

17

Fetita, C., Y.-W. Kim, S. R. Tarando, and P.-Y. Brillet (2019). "Image Biomarkers for Quantitative Analysis of Idiopathic Interstitial Pneumonia". In: *Medical Imaging 2019: Computer-Aided Diagnosis*. Ed. by H. K. Hahn and K. Mori. San Diego, United States: SPIE, p. 44. DOI: 10.1117/12.2511847.

18

19

20

21

Friswell, M. I., J. E. Mottershead, and H. Ahmadian (2001). "Finite-Element Model Updating Using Experimental Test Data: Parametrization and Regularization". In: *Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences* 359.1778. Ed. by N. A. J. Lieven and D. J. Ewins, pp. 169-186. DOI: 10.1098/rsta.2000.0719.

22

23

24

25

Genet, M. (2019). "A Relaxed Growth Modeling Framework for Controlling Growth-Induced Residual Stresses". In: *Clinical Biomechanics* 70, pp. 270-277. DOI: 10.1016/j.clinbiomech.2019.08.015.

26

27

28

Genet, M., L. C. Lee, B. Baillargeon, J. M. Guccione, and E. Kuhl (2016). "Modeling Pathologies of Diastolic and Systolic Heart Failure". In: *Annals of Biomedical Engineering* 44.1, pp. 112-127. DOI: 10.1007/s10439-015-1351-2.

29

30

31

Genet, M., M. Rausch, L. Lee, S. Choy, X. Zhao, G. Kassab, S. Kozerke, J. Guccione, and E. Kuhl (2015a). "Heterogeneous Growth-Induced Prestrain in the Heart". In: *Journal of Biomechanics* 48.10, pp. 2080-2089. DOI: 10.1016/j.jbiomech.2015.03.012.

32

33

34

Genet, M., C. Stoeck, C. Von Deuster, L. Lee, and S. Kozerke (2018). "Equilibrated Warping: Finite Element Image Registration with Finite Strain Equilibrium Gap Regularization". In: *Medical Image Analysis* 50, pp. 1-22. DOI: 10.1016/j.media.2018.07.007.

35

36

37

Genet, M. (2023). "Finite Strain Formulation of the Discrete Equilibrium Gap Principle: Application to Mechanically Consistent Regularization for Large Motion Tracking". In: *Comptes Rendus. Mécanique* 351.G2, pp. 429-458. DOI: 10.5802/crmeca.228.

38

39

40

Genet, M., L. C. Lee, R. Nguyen, H. Haraldsson, G. Acevedo-Bolton, Z. Zhang, L. Ge, K. Ordovas, S. Kozerke, and J. M. Guccione (2014). "Distribution of Normal Human Left Ventricular Myofiber Stress at End Diastole and End Systole: A Target for in Silico Design of Heart Failure Treatments". In: *Journal of Applied Physiology* 117.2, pp. 142-152. DOI: 10.1152/japplphysio.1.00255.2014.

41

42

43

44

45

Genet, M., C. L. Lik, L. Ge, G. Acevedo-Bolton, N. Jeung, A. Martin, N. Cambronero, A. Boyle, Y. Yeghiazarians, and S. Kozerke (2015b). "A Novel Method for Quantifying Smooth Regional Variations in Myocardial Contractility within an Infarcted Human Left Ventricle Based on Delay-Enhanced Magnetic Resonance Imaging". In: *Journal of biomechanical engineering* 137.8, p. 081009.

46

47

48

Genet, M., C. Patte, M. Manoochehrtayebi, and A. Peyraut (2024). *Dolfin_mech*. Zenodo. DOI: 10.5281/ZENODO.10533833.

51

52

Ghattas, O. and K. Willcox (2021). "Learning Physics-Based Models from Data: Perspectives from Inverse Problems and Model Reduction". In: *Acta Numerica* 30, pp. 445-554. DOI: 10.1017/S0962492921000064.

53

54
1

Gibson, G. and N. Pride (1976). "Lung Distensibility. The Static Pressure-Volume Curve of the Lungs and Its Use in Clinical Assessment". In: *British Journal of Diseases of the Chest* 70, pp. 143-184. DOI: 10.1016/0007-0971(76)90027-9.

2

3

Gonsard, A., M. Genet, and D. Drummond (2024). "Digital Twins for Chronic Lung Diseases". In: *European Respiratory Review* 33.174, p. 240159. DOI: 10.1183/16000617.0159-2024.

4

7

Grédiac, M., F. Pierron, S. Avril, and E. Toussaint (2006). "The Virtual Fields Method for Extracting Constitutive Parameters From Full-Field Measurements: A Review". In: *Strain* 42.4, pp. 233-253. DOI: 10.1111/j.1475-1305.2006.tb01504.x.

8

Grédiac, M., F. Hild, and A. Pineau, eds. (2013). *Full-Field Measurements and Identification in Solid Mechanics*. 1st ed. Wiley. DOI: 10.1002/9781118578469.

10

Hall, J. K., J. H. T. Bates, R. Krishnan, J. H. Kim, Y. Deng, K. R. Lutchen, and B. Suki (2024). "Elucidating the Interaction between Stretch and Stiffness Using an Agent-Based Spring Network Model of Progressive Pulmonary Fibrosis". In: *Frontiers in Network Physiology* 4, p. 1396383. DOI: 10.3389/fnetp.2024.1396383.

12

13

14

Hinz, B. and B. Suki (2016). "Does Breathing Amplify Fibrosis?" In: *American Journal of Respiratory and Critical Care Medicine* 194.1, pp. 9-11. DOI: 10.1164/rccm.201601-0149ED.

16

Hofmann, W., B. Asgharian, and R. Winkler-Heil (2002). "Modeling Intersubject Variability of Particle Deposition in Human Lungs". In: *Journal of Aerosol Science* 33.2, pp. 219-235. DOI: 10.1016/S0021-8502(01)00167-7.

18

19

Hu, S., E. Hoffman, and J. Reinhardt (2001). "Automatic Lung Segmentation for Accurate Quantitation of Volumetric X-ray CT Images". In: *IEEE Transactions on Medical Imaging* 20.6, pp. 490-498. DOI: 10.1109/42.929615.

21

22

Hua, X., Y. Ni, and J. Ko (2009). "Adaptive Regularization Parameter Optimization in Output-Error-Based Finite Element Model Updating". In: *Mechanical Systems and Signal Processing* 23.3, pp. 563-579. DOI: 10.1016/j.ymssp.2008.05.002.

24

25

Hurtado, D. E., N. Villarroel, C. Andrade, J. Retamal, G. Bugedo, and A. Bruhn (2017). "Spatial Patterns and Frequency Distributions of Regional Deformation in the Healthy Human Lung". In: *Biomechanics and modeling in mechanobiology* 16, pp. 1413-1423.

27

28

Kallet, R. H. (2015). "A Comprehensive Review of Prone Position in ARDS". In: *Respiratory care* 60.11, pp. 1660-1687.

30

Kaneko, K., J. Milic-Emili, M. B. Dolovich, A. Dawson, and D. V. Bates (1966). "Regional Distribution of Ventilation and Perfusion as a Function of Body Position." In: *Journal of Applied Physiology* 21.3, pp. 767-777.

32

33

Kaul, H. (2019). "Respiratory Healthcare by Design: Computational Approaches Bringing Respiratory Precision and Personalised Medicine Closer to Bedside". In: *Morphologie* 103.343, pp. 194-202.

36

Kavanagh, K. T. and R. W. Clough (1971). "Finite Element Applications in the Characterization of Elastic Solids". In: *International Journal of Solids and Structures* 7.1, pp. 11-23. DOI: 10.1016/0020-7683(71)90015-1.

38

39

Kerfriden, P., A. Kundu, and S. Claus (2019). "Adaptivity in Bayesian Inverse Finite Element Problems: Learning and Simultaneous Control of Discretisation and Sampling Errors". In: *Materials* 12.4, p. 642. DOI: 10.3390/ma12040642.

42

43

Laville, C., C. Fetita, T. Gille, P.-Y. Brillet, H. Nunes, J.-F. Bernaudin, and M. Genet (2023). "Comparison of Optimization Parametrizations for Regional Lung Compliance Estimation Using Personalized Pulmonary Poromechanical Modeling". In: *Biomechanics and Modeling in Mechanobiology*, pp. 1-13.

45

46

Lazarus, A., D. Dalton, D. Husmeier, and H. Gao (2022). "Sensitivity Analysis and Inverse Uncertainty Quantification for the Left Ventricular Passive Mechanics". In: *Biomechanics and Modeling in Mechanobiology* 21.3, pp. 953-982. DOI: 10.1007/s10237-022-01571-8.

48

49

Le Gratiet, L., C. Cannamela, and B. Iooss (2014). "A Bayesian Approach for Global Sensitivity Analysis of (Multifidelity) Computer Codes". In: *SIAM/ASA Journal on Uncertainty Quantification* 2.1, pp. 336-363. DOI: 10.1137/130926869.

51

52

Leclerc, H., J.-N. Périé, S. Roux, and F. Hild (2011). "Voxel-Scale Digital Volume Correlation". In: *Experimental Mechanics* 51.4, pp. 479-490. DOI: 10.1007/s11340-010-9407-6.

55

56

Lee, L. C. and M. Genet (2019). "Validation of Equilibrated Warping—Image Registration with Mechanical Regularization—On 3D Ultrasound Images". In: *Functional Imaging and Modeling*
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

of the Heart. Ed. by Y. Coudière, V. Ozenne, E. Vigmond, and N. Zemzemi. Vol. 11504. Cham: Springer International Publishing, pp. 334–341. DOI: 10.1007/978-3-030-21949-9_36.

Lee, L. C., M. Genet, A. B. Dang, L. Ge, J. M. Guccione, and M. B. Ratcliffe (2014). “Applications of Computational Modeling in Cardiac Surgery”. In: *Journal of Cardiac Surgery: Including Mechanical and Biological Support for the Heart and Lungs* 29.3, pp. 293–302.

Li, L., H. Li, J. Shu, J. Pan, X. Chen, M. Ying, Y. Xu, D. Wang, and P. Pang (2019). “Influences of Scan-Position on Clinical Ultra-High-Resolution CT Scanning: A Preliminary Study”. In: *Scientific Reports* 9.1, p. 1134. DOI: 10.1038/s41598-018-37514-6.

Logg, A., K.-A. Mardal, and G. Wells, eds. (2012). *Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book*. Vol. 84. Lecture Notes in Computational Science and Engineering. Berlin, Heidelberg: Springer Berlin Heidelberg. DOI: 10.1007/978-3-642-23099-8.

Maghsoudi-Ganjeh, M., C. A. Mariano, S. Sattari, H. Arora, and M. Eskandari (2021). “Developing a Lung Model in the Age of COVID-19: A Digital Image Correlation and Inverse Finite Element Analysis Framework”. In: *Frontiers in Bioengineering and Biotechnology* 9, p. 684778. DOI: 10.3389/fbioe.2021.684778.

Marquis, A. D., A. Arnold, C. Dean-Bernhoft, B. E. Carlson, and M. S. Olufsen (2018). “Practical Identifiability and Uncertainty Quantification of a Pulsatile Cardiovascular Model”. In: *Mathematical Biosciences* 304, pp. 9–24. DOI: 10.1016/j.mbs.2018.07.001.

Mendoza, A., J. Neggers, F. Hild, and S. Roux (2019). “Complete Mechanical Regularization Applied to Digital Image and Volume Correlation”. In: *Computer Methods in Applied Mechanics and Engineering* 355, pp. 27–43. DOI: 10.1016/j.cma.2019.06.005.

Milic-Emili, J., J. Mead, J. M. Turner, and E. M. Glauser (1964). “Improved Technique for Estimating Pleural Pressure from Esophageal Balloons”. In: *Journal of Applied Physiology* 19.2, pp. 207–211. DOI: 10.1152/jappl.1964.19.2.207.

Miller, M. I., G. E. Christensen, Y. Amit, and U. Grenander (1993). “Mathematical Textbook of Deformable Neuroanatomies.” In: *Proceedings of the National Academy of Sciences* 90.24, pp. 11944–11948. DOI: 10.1073/pnas.90.24.11944.

Mraz, T., S. Asgari, A. Karimi, M.-K. Breyer, S. Hartl, O. Sunanta, A. Ofenheimer, O. C. Burghuber, A. Zacharasiewicz, B. Lamprecht, C. Schiffers, E. F. M. Wouters, and R. Breyer-Kohansal (2024). “Updated Reference Values for Static Lung Volumes from a Healthy Population in Austria”. In: *Respiratory Research* 25.1, p. 155. DOI: 10.1186/s12931-024-02782-6.

Murata, N., S. Ito, K. Furuya, N. Takahara, K. Naruse, H. Aso, M. Kondo, M. Sokabe, and Y. Hasegawa (2014). “Ca2+ Influx and ATP Release Mediated by Mechanical Stretch in Human Lung Fibroblasts”. In: *Biochemical and Biophysical Research Communications* 453.1, pp. 101–105. DOI: 10.1016/j.bbrc.2014.09.063.

Myers, C. A., P. J. Laz, K. B. Shelburne, and B. S. Davidson (2015). “A Probabilistic Approach to Quantify the Impact of Uncertainty Propagation in Musculoskeletal Simulations”. In: *Annals of Biomedical Engineering* 43.5, pp. 1098–1111. DOI: 10.1007/s10439-014-1181-7.

Nagel, J. B. and B. Sudret (2016). “A Unified Framework for Multilevel Uncertainty Quantification in Bayesian Inverse Problems”. In: *Probabilistic Engineering Mechanics* 43, pp. 68–84. DOI: 10.1016/j.probengmech.2015.09.007.

Patte, C. (2020). “Personalized Pulmonary Mechanics: Modeling, Estimation and Application to Pulmonary Fibrosis”. PhD thesis. Institut polytechnique de Paris.

Patte, C., P.-Y. Brillet, C. Fetita, J.-F. Bernaudin, T. Gille, H. Nunes, D. Chapelle, and M. Genet (2022a). “Estimation of Regional Pulmonary Compliance in Idiopathic Pulmonary Fibrosis Based on Personalized Lung Poromechanical Modeling”. In: *Journal of Biomechanical Engineering* 144.9, p. 091008. DOI: 10.1115/1.4054106.

Patte, C., M. Genet, and D. Chapelle (2022b). “A Quasi-Static Poromechanical Model of the Lungs”. In: *Biomechanics and Modeling in Mechanobiology* 21.2, pp. 527–551. DOI: 10.1007/s10237-021-01547-0.

Pellegrino, D., S. Casas-Recasens, R. Faner, P. Palange, and A. Agusti (2023). “When GETomics Meets Aging and Exercise in COPD”. In: *Respiratory Medicine* 216, p. 107294. DOI: 10.1016/j.rmed.2023.107294.
Peyraut, A. and M. Genet (2024a). "A Model of Mechanical Loading of the Lungs Including Gravity and a Balancing Heterogeneous Pleural Pressure". In: *Biomechanics and Modeling in Mechanobiology*. DOI: 10.1007/s10237-024-01876-w.

— (2024b). "A Model of Mechanical Loading of the Lungs Including Gravity and a Balancing Heterogeneous Pleural Pressure". In.

— (2024c). *Uncertainty Quantification Paper Demos*. Zenodo. DOI: 10.5281/ZENODO.13928262.

— (2025). "Finite Strain Formulation of the Discrete Equilibrium Gap Principle: Application to Direct Parameter Estimation from Large Full-Fields Measurements". In: *Comptes Rendus. Mécanique* 353.G1, pp. 259-308. DOI: 10.5802/crmeca.279.

Plantier, L., A. Cazes, A.-T. Dinh-Xuan, C. Bancal, S. Marchand-Adam, and B. Crestani (2018). "Physiology of the Lung in Idiopathic Pulmonary Fibrosis". In: *European Respiratory Review* 27.147, p. 170062. DOI: 10.1183/16000617.0062-2017.

Rausch, M. K., A. M. Zöllner, M. Genet, B. Baillargeon, W. Bothe, and E. Kuhl (2017). "A Virtual Sizing Tool for Mitral Valve Annuloplasty". In: *International Journal for Numerical Methods in Biomedical Engineering* 33.2, e02788. DOI: 10.1002/cnm.2788.

Rossi, M., P. Lava, F. Pierron, D. Debruyne, and M. Sasso (2015). "Effect of DIC Spatial Resolution, Noise and Interpolation Error on Identification Results with the VFM". In: *Strain* 51.3, pp. 206-222. DOI: 10.1111/str.12134.

Roth, C. J., T. Becher, I. Frerichs, N. Weiler, and W. A. Wall (2017). "Coupling of EIT with Computational Lung Modeling for Predicting Patient-Specific Ventilatory Responses". In: *Journal of Applied Physiology* 122.4, pp. 855-867. DOI: 10.1152/japplphysiol.00236.2016.

Roux, S. and F. Hild (2020). "Optimal Procedure for the Identification of Constitutive Parameters from Experimentally Measured Displacement Fields". In: *International Journal of Solids and Structures* 184, pp. 14-23. DOI: 10.1016/j.ijsolstr.2018.11.008.

Saltelli, A. and I. M. Sobol (1995). "Sensitivity Analysis for Nonlinear Mathematical Models: Numerical Experience". In: *Matematicheskoe Modelirovanie* 7.11, pp. 16-28.

Santos, V. J., C. D. Bustamante, and F. J. Valero-Cuevas (2009). "Improving the Fitness of High-Dimensional Biomechanical Models via Data-Driven Stochastic Exploration". In: *IEEE Transactions on Biomedical Engineering* 56.3, pp. 552-564. DOI: 10.1109/TBME.2008.2006033.

Sicard, D., A. J. Haak, K. M. Choi, A. R. Craig, L. E. Fredenburgh, and D. J. Tschumperlin (2018). "Aging and Anatomical Variations in Lung Tissue Stiffness". In: *American Journal of Physiology-Lung Cellular and Molecular Physiology* 314.6, pp. L946-L955. DOI: 10.1152/ajplung.00415.2017.

Siekmann, I., J. Sneyd, and E. J. Crampin (2012). "MCMC Can Detect Nonidentifiable Models". In: *Biophysical Journal* 103.11, pp. 2275-2286. DOI: 10.1016/j.bpj.2012.10.024.

Škardová, K., M. Rambausek, R. Chabiniok, and M. Genet (2019). "Mechanical and Imaging Models-Based Image Registration". In: *Lecture Notes in Computational Vision and Biomechanics*. Cham: Springer International Publishing, pp. 77-85. DOI: 10.1007/978-3-030-32040-9_9.

Smith, N., A. de Vecchi, M. McCormick, D. Nordsletten, O. Camara, A. F. Frangi, H. Delingette, M. Sermesant, J. Relan, and N. Ayache (2011). "euHeart: Personalized and Integrated Cardiac Care Using Patient-Specific Cardiovascular Modelling". In: *Interface focus* 1.3, pp. 349-364.

Sobol, I., S. Tarantola, D. Gatelli, S. Kucherenko, and W. Mauntz (2007). "Estimating the Approximation Error When Fixing Unessential Factors in Global Sensitivity Analysis". In: *Reliability Engineering & System Safety* 92.7, pp. 957-960. DOI: 10.1016/j.ress.2006.07.001.

Sutton, M. A. and F. Hild (2015). "Recent Advances and Perspectives in Digital Image Correlation". In: *Experimental Mechanics* 55.1, pp. 1-8. DOI: 10.1007/s11340-015-9991-6.

Tarantola, A. (2005). *Inverse Problem Theory and Methods for Model Parameter Estimation*. SIAM.

Tawhai, M. H., M. P. Nash, C.-L. Lin, and E. A. Hoffman (2009). "Supine and Prone Differences in Regional Lung Density and Pleural Pressure Gradients in the Human Lung with Constant Shape". In: *Journal of Applied Physiology* 107.3, pp. 912-920.

Thannickal, V. J., C. A. Henke, J. C. Horowitz, P. W. Noble, J. Roman, P. J. Sime, Y. Zhou, R. G. Wells, E. S. White, and D. J. Tschumperlin (2014). "Matrix Biology of Idiopathic Pulmonary Fibrosis". In: *The American Journal of Pathology* 184.6, pp. 1643-1651. DOI: 10.1016/j.ajpat.h.2014.02.003.

Tosin, M., A. M. A. Côrtes, and A. Cunha (2020). "A Tutorial on Sobol' Global Sensitivity Analysis Applied to Biological Models". In: *Networks in Systems Biology*. Ed. by F. A. B. Da Silva,
1

N. Carels, M. Trindade Dos Santos, and F. J. P. Lopes. Vol. 32. Cham: Springer International Publishing, pp. 93–118. DOI: 10.1007/978-3-030-51862-2_6.

2

3 Veress, A. I., G. T. Gullberg, and J. A. Weiss (2005). “Measurement of Strain in the Left Ventricle during Diastole with Cine-MRI and Deformable Image Registration”. In: *Journal of Biomechanical Engineering* 127.7, pp. 1195–1207. DOI: 10.1115/1.2073677.

4

5

6 West, J. B. and F. L. Matthews (1972). “Stresses, Strains, and Surface Pressures in the Lung Caused by Its Weight”. In: *Journal of Applied Physiology* 32.3, pp. 332–345.

7

8 Wieland, F.-G., A. L. Hauber, M. Rosenblatt, C. Tönsing, and J. Timmer (2021). “On Structural and Practical Identifiability”. In: *Current Opinion in Systems Biology* 25, pp. 60–69. DOI: 10.1016/j.coisb.2021.03.005.

10

11 Wijsenbeek, M. and V. Cottin (2020). “Spectrum of Fibrotic Lung Diseases”. In: *New England Journal of Medicine* 383.10, pp. 958–968.

12

13 Wirthl, B., S. Brandstaeter, J. Nitzler, B. A. Schrefler, and W. A. Wall (2023). “Global Sensitivity Analysis Based on Gaussian-process Metamodelling for Complex Biomechanical Problems”. In: *International Journal for Numerical Methods in Biomedical Engineering* 39.3, e3675. DOI: 10.1002/cnm.3675.

14

15

16

17 Wu, H., Y. Yu, H. Huang, Y. Hu, S. Fu, Z. Wang, M. Shi, X. Zhao, J. Yuan, J. Li, et al. (2020). “Progressive Pulmonary Fibrosis Is Caused by Elevated Mechanical Tension on Alveolar Stem Cells”. In: *Cell* 180.1, pp. 107–121.

19

20 Wu, X., T. Kozlowski, H. Meidani, and K. Shirvan (2018). “Inverse Uncertainty Quantification Using the Modular Bayesian Approach Based on Gaussian Process, Part 1: Theory”. In: *Nuclear Engineering and Design* 335, pp. 339–355. DOI: 10.1016/j.nucengdes.2018.06.004.

21

22

23 Wu, X., K. Shirvan, and T. Kozlowski (2019). “Demonstration of the Relationship between Sensitivity and Identifiability for Inverse Uncertainty Quantification”. In: *Journal of Computational Physics* 396, pp. 12–30. DOI: 10.1016/j.jcp.2019.06.032.

24

25

26 Ying, Z. (1991). “Asymptotic Properties of a Maximum Likelihood Estimator with Data from a Gaussian Process”. In: *Journal of Multivariate Analysis* 36.2, pp. 280–296. DOI: 10.1016/0047-259X(91)90062-7.

27

28

29 Zhang, H. (2004). “Inconsistent Estimation and Asymptotically Equal Interpolations in Model-Based Geostatistics”. In: *Journal of the American Statistical Association* 99.465, pp. 250–261. DOI: 10.1198/0162145040000000241.

30

31
# A Recall of the equations of the model

## A.1 Poromechanical framework

The chosen working framework is based on the clinical data we have access to. As lungs are a porous medium, and as a porosity field can be extracted from the CT-scans, we choose to work in a poromechanical framework. Additionally, our quasi-static problem can be decomposed into two computations: the unloading problem and the loading problem, for which kinematics operators are introduced hereafter.

**Kinematics of the loading problem** The loading problem corresponds to the computation of the end-inhalation configuration from the unloaded configuration. In this article, quantities written in capital letters are associated to the unloaded configuration, while lower case letters refer to quantities associated to a deformed configuration. Let us first define quantities and operators associated to the loading problem.

The function $\chi$ represents the mapping from a position $\underline{X}$ of a material point in the unloaded configuration, denoted $\Omega$, to the position $\underline{x}$ of the same material point in the deformed configuration, denoted $\omega$, i.e.,

$$ \underline{\chi} : \underline{X} \mapsto \underline{x}. \qquad (11) $$

The mapping function is involved in the definition of the displacement field from the unloaded to the deformed configuration $\underline{U}$,

$$ \underline{U} := \underline{\chi}(\underline{X}) - \underline{X} = \underline{x} - \underline{X}, \qquad (12) $$

as well as in the definition of the deformation gradient tensor $\underline{F}$,

$$ \underline{F} := \underline{\nabla} \underline{\chi} = \underline{1} + \underline{\nabla} \underline{U}. \qquad (13) $$

Different quantities derive directly from the deformation gradient tensor, such as the volume change $J$, defined as

$$ J := \det(\underline{F}), \qquad (14) $$

the right Cauchy tensor, expressed as

$$ \underline{C} := \underline{F}^T \cdot \underline{F}, \qquad (15) $$

and the Green-Lagrange strain tensor, defined as

$$ \underline{E} := \frac{1}{2}(\underline{C} - \underline{1}) = \frac{1}{2}(\underline{F}^T \cdot \underline{F} - \underline{1}). \qquad (16) $$

We also call $\underline{N}$ the outer normal in the unloaded configuration, and $\underline{\Gamma} := \partial\Omega$ the unloaded configuration boundary, linked to the outer normal in the deformed configuration $\underline{n}$ through the following relationship:

$$ \underline{n} = \frac{\underline{F}^{-T} \underline{N}}{\|\underline{F}^{-T} \underline{N}\|}. \qquad (17) $$

Note that the loaded configuration boundary is denoted $\underline{\gamma} := \partial\omega$

**Kinematics of the unloading problem** As for the loading problem, different operators can be defined to describe the kinematics of the unloading problem, i.e., the computation of the unknown unloaded configuration from a known deformed configuration.

The mapping function $\underline{\chi}$ defined in Equation (11) allows us to define the displacement field from the deformed configuration to the unloaded configuration, written $\underline{u}$:

$$ \underline{u}(\underline{x}) := \underline{\chi}^{-1}(\underline{x}) - \underline{x} = \underline{X} - \underline{x}, \qquad (18) $$
as well as the deformation gradient tensor from the deformed configuration to the unloaded configuration $\underline{f}$:

$$ \underline{f} := \underline{\nabla} \chi^{-1}. \qquad (19) $$

Similarly to the loading problem, the deformation gradient tensor $\underline{f}$ can also be used to define the volume change $j$ between the deformed configuration and the unloaded configuration,

$$ j := \det(\underline{f}). \qquad (20) $$

**Poromechanics operators** As a reminder, we are working in a poromechanical framework, which means that each point of the domain contains both a fluid –the air– and a solid –tissues and blood– phase. Quantities associated to the fluid phase are written with the subscript “f” and quantities associated to the solid phase are written with the subscript “s”.

Let us consider an infinitesimal volume of matter $|d\Omega|$ in the unloaded configuration $\Omega$. $|d\Omega|$ is composed of an infinitesimal volume of fluid $|d\Omega_f|$, and of an infinitesimal volume of solid $|d\Omega_s|$, such that

$$ |d\Omega| = |d\Omega_f| + |d\Omega_s|. \qquad (21) $$

In order to describe the fluid and solid volume fraction in the unloaded configuration, we can therefore introduce two variables, $\Phi_{f0}$ (called the reference porosity) and $\Phi_{s0}$, defined as:

$$ \left\{ \begin{aligned} \Phi_{f0} &:= \frac{|d\Omega_{f0}|}{|d\Omega|} \\ \Phi_{s0} &:= \frac{|d\Omega_{s0}|}{|d\Omega|} \end{aligned} \right. . \qquad (22) $$

Similarly, in the deformed configuration, an infinitesimal volume $|d\omega|$ is composed of an infinitesimal volume of fluid $|d\omega_f|$, and of an infinitesimal volume of solid $|d\omega_s|$, such that

$$ |d\omega| = |d\omega_f| + |d\omega_s|. \qquad (23) $$

We can therefore introduce two variables $\phi_f$ (called the deformed porosity) and $\phi_s$, describing respectively the fluid and solid proportion in the deformed configuration. These variables are defined as:

$$ \left\{ \begin{aligned} \phi_f &:= \frac{|d\omega_f|}{|d\omega|} \\ \phi_s &:= \frac{|d\omega_s|}{|d\omega|} \end{aligned} \right. . \qquad (24) $$

We can also introduce additional constituent volume fractions, $\phi_{f0}$ and $\phi_{s0}$, linking the unloaded fluid and solid fractions and the deformed volume, as well as constituent volume fractions, $\Phi_f$ and $\Phi_s$ linking the deformed fluid and solid fractions and the unloaded volume:

$$ \left\{ \begin{aligned} \phi_{f0} &:= \frac{|d\Omega_{f0}|}{|d\omega|} \\ \phi_{s0} &:= \frac{|d\Omega_{s0}|}{|d\omega|} \end{aligned} \right. \quad \text{and} \quad \left\{ \begin{aligned} \Phi_f &:= \frac{|d\omega_f|}{|d\Omega|} \\ \Phi_s &:= \frac{|d\omega_s|}{|d\Omega|} \end{aligned} \right. . \qquad (25) $$

## A.2 Equilibrium equations

**Stresses** Let us first introduce two usual stress tensors, used for writing the equilibrium equations: $\underline{\sigma}$, the Cauchy stress tensor, defined in the deformed configuration, and $\underline{\Sigma}$, the second Piola-Kirchhoff stress tensor, defined in the reference configuration. The second Piola-Kirchhoff stress tensor is linked to the Cauchy stress tensor through the following relation (the mapping between deformed and unloaded configurations does not appear for the sake of simplicity):

$$ \underline{\Sigma} := J \underline{f}^{-1} \underline{\sigma} \underline{f}^{-T}, \qquad (26) $$

Conversely, the Cauchy stress tensor can be expressed in relation to the second Piola-Kirchhoff stress tensor:

$$ \underline{\sigma} = j \underline{f}^{-1} \underline{\Sigma} \underline{f}^{-T}. \qquad (27) $$
1 **Equilibrium** We express the equilibrium via the principle of virtual work, which derives directly
2 from the principle of conservation of momentum:

$$
\underline{\underline{\sigma}} / \underline{\underline{\Sigma}} \mid W_{\text{int}}(\underline{\underline{\sigma}} / \underline{\underline{\Sigma}} ; \underline{\underline{U}}^{*}) = W_{\text{ext}}(\underline{\underline{\sigma}} / \underline{\underline{\Sigma}} ; \underline{\underline{U}}^{*}) \forall \underline{\underline{U}}^{*}, \quad (28)
$$

3 where U* represents any kinematically admissible to zero virtual field. The expressions of the
4 internal and external virtual work for the loading and the unloading problems are detailed hereafter.

5 **Virtual work of internal forces** For the unloading problem, the virtual work of internal forces
6 can be expressed as

$$
W_{\text{int}}(\underline{u}; \underline{u}^*) := \int_{\omega} \underline{\underline{\sigma}}(\underline{u}) : \underline{\underline{\epsilon}}(\underline{u}^*) d\omega, \qquad (29)
$$

7 , where $\underline{\underline{\epsilon}}$ represents the strain tensor, defined as

$$
\underline{\underline{\epsilon}} := \frac{1}{2} (\underline{\underline{\nabla}} \underline{u} + \underline{\underline{\nabla}} \underline{u}^{\mathrm{T}}). \qquad (30)
$$

9 In the case of the loading problem, it is expressed as

$$
W_{\text{int}}(\underline{U}; \underline{U}^*) := \int_{\Omega} \underline{\underline{\Sigma}}(\underline{U}) : d\underline{U} \underline{\underline{E}} \cdot \underline{U}^* d\Omega. \quad (31)
$$

10 where dU E · U* is the first differential of the Green-Lagrange strain tensor. Note that this equation
11 includes implicitly the internal fluid pressure field, i.e., the forces applied by airway structures to
12 the parenchyma, as detailed in [Patte et al. 2022b].

13 **Virtual work of external forces** Two main loads, gravity -denoted as g- and the pleural
14 pressure p̃_pl, have an impact on the lung behavior. However, if the pleural pressure is taken
15 into account in most macroscale models, it is not the case of gravity [Patte et al. 2022b]. The
16 effects of gravity are indeed small compared to those of the pleural pressure. Additionally, the
17 implementation of gravity into a model supposes the formulation of physiological counterbalancing
18 boundary conditions, which is not straightforward.

19 However, there are many indicators that gravity does have an impact on the physiology of the
20 lungs [West et al. 1972; Tawhai et al. 2009; Kallet 2015]. In order to implement gravity into the
21 model, physiological boundary conditions have to be formulated. We can consider that gravity is
22 counterbalanced in the lungs by the pleural pressure and by reaction forces of the surroundings. If
23 modeling explicitely the surroundings is, as for many biomechanical problems [Chabiniok et al. 2012],
24 too complex, we can approximate all these reaction forces as an heterogeneous counterbalancing
25 pressure p<sub>pl</sub> applied on the outer lung surface. In order to be physiological, this pressure needs to be
26 as regular as possible. As this is a counterbalancing pressure, it should aslo verify the equilibrium
27 of forces and momentum.

28 As there already exists a physiological pressure p̃_pl, which corresponds to the subatmospheric
29 pressure in the pleural cavity, we choose to formulate p<sub>pl</sub> as a counterbalancing pressure, which
30 also includes the pressure p̃_pl. p<sub>pl</sub> can therefore be obtained through the following constrained
31 minimization, as detailed in [Peyraut et al. 2024a]:

$$
p_{\text{pl}} = \begin{array}{ll} \underset{\text{p such that}}{\operatorname{argmin}} & \left\{ \frac{1}{2} \int_{\gamma} \left( \underline{\underline{\nabla}}_{\text{s}} p - \underline{\underline{\nabla}}_{\text{s}} \tilde{p}_{\text{pl}} \right)^2 d\gamma \right\}, \quad (32) \\ & \left\{ \begin{array}{l} \int_{\omega} \rho \phi_{\text{s}} \underline{\underline{g}} \, d\omega - \int_{\gamma} \underline{\underline{p}}_{\text{n}} \, d\gamma = \underline{\underline{0}} \\ \int_{\gamma} \tilde{\underline{\underline{x}}} \times \underline{\underline{p}}_{\text{n}} \, d\gamma = \underline{\underline{0}} \\ \int_{\gamma} (\underline{\underline{p}} - \tilde{\underline{\underline{p}}}_{\text{pl}}) \, d\gamma = \underline{\underline{0}} \end{array} \right\} \end{array}
$$

32 where $\tilde{\underline{x}} := \begin{pmatrix} \tilde{x} \\ \tilde{y} \\ \tilde{z} \end{pmatrix}$ corresponds to the center of mass of the structure. For the implementation in our
33 model, $p_{\text{pl}}$ is obtained using Lagrange multipliers. Details are provided in [Peyraut et al. 2024a],
$$ \tilde{p}_{pl} = \tilde{p}_0 + K(z - \bar{z}), \quad (33) $$

$$ W_{\text{ext}}(\underline{u}; \underline{u}^*) = \int_{\omega} \rho \phi_s \underline{g} \cdot \underline{u}^* d\omega + \int_{\gamma} p_{\text{pl}} \underline{n} \cdot \underline{u}^* d\gamma, \quad (34) $$

$$ W_{\text{ext}}(\underline{U}; \underline{U}^*) = \int_{\Omega} \rho \Phi_{s0} \underline{g} \cdot \underline{U}^* d\Omega + \int_{\Gamma} p_{\text{pl}} \underline{J} \underline{F}^{-T} \cdot \underline{N} \cdot \underline{U}^* d\Gamma. \quad (35) $$

## A.3 Constitutive framework

$$ \Psi = \Psi_s + \Psi_f. \quad (36) $$

$$ \underline{\Sigma} = \frac{\partial \Psi}{\partial \underline{E}} = \underline{\Sigma}_s - p_f \underline{J} \underline{C}^{-1}, \quad (37) $$

$$ \underline{\Sigma}_s := \frac{\partial \Psi_s}{\partial \underline{E}}, \quad (38) $$

$$ p_f = -\frac{\partial \Psi_s}{\partial \Phi_s}. \quad (39) $$

$$ \Psi_s = \Psi_{\text{skel}} + \Psi_{\text{bulk}} + \Psi_{\text{por}}. \quad (40) $$

$$ \Psi_{\text{skel}} = \alpha(e^{\delta(J^2 - 2 \ln J)} - 1) + \beta_1(I_1 - 3 - 2 \ln(J)) + \beta_2(I_2 - 3 - 4 \ln(J)), \quad (41) $$

<footnote n="1">1 and the formulation is also detailed in Appendix C. Additionally, following the observations of
2 [Kaneko et al. 1966; Barrau 2024] suggesting that the pleural pressure has the form of a gradient,
3 we chose to express, the pleural pressure $\tilde{p}_{pl}$ as:</footnote>

<footnote n="4">4 where $\tilde{p}_0$ and K are constants, with $K = \rho_{plg}$, $\rho_{pl}$ being the density of the pleural liquid and g the acceleration of gravity, which is a positive constant for a patient in supine position and a negative constant for a patient in prone position.</footnote>

<footnote n="7">7 We can therefore express the virtual work of external forces of the unloading problem as</footnote>

<footnote n="8">8 and the virtual work of external forces of the loading problem as</footnote>

<footnote n="10">10 **Behavior** Let us call $\Psi$ the system Helmholtz free energy. Following [Coussy 2004; Chapelle et al. 2014], we split $\Psi$ into two contributions: a solid contribution, $\Psi_s$, and a fluid contribution $\Psi_f$, such that</footnote>

<footnote n="13">13 Note that $\Psi_s = \Psi_s(E, \Phi_f)$ and $\Psi_f = \Psi_f(\Phi_f)$, and that $\Phi_f$ does not depend on $E$. Detailed explanations can be found [Patte 2020]. Following the second principle of thermodynamics, we can then link the second Piola-Kirchhoff stress tensor to the system Helmholtz free energy $\Psi$ through the relation:</footnote>

<footnote n="17">17 where $\underline{\Sigma}_s$ is the solid stress defined as</footnote>

<footnote n="18">18 and $p_f$ is the internal fluid pressure field, expressed as</footnote>

<footnote n="19">19 This last equation represents the internal equilibrium between the solid and the fluid phase. The second Piola-Kirchhoff stress tensor therefore only depends on $\Psi_s$ –i.e., no dependency on $\Psi_f$–. As a result, only the solid free energy needs to be formulated.</footnote>

<footnote n="22">22 **Solid free energy** As for the system free energy, the solid free energy can be split into three contributions [Coussy 2004; Chapelle et al. 2014]. The first contribution, $\Psi_{skel}$ represents the behavior of the solid phase, while the second contribution $\Psi_{bulk}$ represents the changes in porosity. Lastly, the potential $\Psi_{por}$ ensures positive –and hence, physiological– values of porosity. The solid free energy can therefore be written as:</footnote>

<footnote n="27">27 The skeleton free energy $\Psi_{skel}$ is based on an exponential, which is classical for soft tissues [Demiray 1972; Caruel et al. 2014; Genet et al. 2015b]:</footnote>
1 where $I_1$ and $I_2$ are the two first invariants of $\mathcal{C}$, and $\alpha$ and $\delta$ are material constants. The bulk free

2 energy $\Psi_{\text{bulk}}$ is then expressed as

$$ \Psi_{\text{bulk}} = \kappa \left( \left( \frac{\Phi_s}{\Phi_{s0}} \right)^2 - 1 - 2 \ln \left( \frac{\Phi_s}{\Phi_{s0}} \right) \right), \quad (42) $$

3 where $\beta_1$ and $\kappa$ are material constants. Lastly, the pore energy $\Psi_{\text{por}}$ is expressed as

$$ \Psi_{\text{por}} = \zeta \begin{cases} +\infty & \text{if } \Phi_f < 0 \text{ or } \Phi_{f0} < 0 \\ (\Phi_{f0}/\Phi_f - 1)^2 & \text{if } \Phi_f \le \Phi_{f0} \\ (\Phi_f/\Phi_{f0} - 1)^2 & \text{if } \Phi_{f0} < \Phi_f \end{cases}, \quad (43) $$

4 where $\zeta$ is a material constant. When the porosity becomes very small, the pore energy increases,

5 preventing the porosity to reach 0, i.e. a full collapse of the pores, or negative values, i.e. non-

6 physiological values.
## Appendix B Using gravity to improve the estimation process: a proof of concept on the linear elastic spring

The idea of using gravity to improve the estimation process can be validated on a very simple and highly identifiable case: a linear spring under tension. Let us consider a simple spring, clamped at one extremity and pulled with a force f at the other. In a first approach, let us neglect the mass of the spring. The displacement of the extremity of the spring is therefore simply given by Hooke's law:

$$u(\alpha, f) = \frac{f}{\alpha} \qquad (44)$$

where α is the stiffness of the spring

Applying the pipeline described in the Methods section is equivalent to consider a measure:

$$u_{\text{meas}} = \frac{f_0}{\alpha_0} + n \qquad (45)$$

where f₀ and α₀ are values of f and α arbitrarily chosen and n is a noise added to u(f₀, α₀), whose value depends on the SNR considered.

If there is no noise, one can expect to be able to identify perfectly the force knowing the stiffness or the stiffness knowing the noise. However, it should be impossible to identify both the force and the stiffness simultaneously. The distribution of the error on the estimated stiffness parameter (assuming the force of tension to be known) and the distribution of the error on the estimated force of tension (assuming the stiffness to be known) are presented in Figure 14 as a function of the Signal to Noise Ratio (SNR):

<figure><img src="image_18.png" /><figcaption>Figure 14: Error distribution on the estimated stiffness knowing the force (left) and on the estimated force knowing the stiffness (right)</figcaption></figure>

In this example, for the figure on the left, the force of tension of the spring is assumed to be known exactly, the displacement is assumed to be measured with a given level of noise, and the estimated parameter is the spring stiffness. For any level of noise, our method provides the estimation error for the stiffness. For large SNRs –i.e., for a small impact of the noise–, the optimization always converges towards the exact parameter value; whereas for small SNRs, the obtained parameters are far away from the ground-truth. The smaller the SNR, the higher the impact of the noise and the higher the error made on the estimated parameter. Our pipeline thus allows us to quantify the identifiability of the stiffness depending on the noise on the measures. These results are consistent with the expected error distribution in the case of a linear spring and therefore provide a first element of validation of our pipeline. The same observation can be made for the figure on the right, for which the estimated parameter is the force of tension of the spring, while the stiffness is assumed to be known exactly.
The distribution of the error on the estimated stiffness parameter and on the estimated force of tension, when both parameters are estimated simultaneously, is presented in Figure 15 as a function of the Signal to Noise Ratio (SNR). For any level of noise, the error distribution on the estimated parameters is wider than the initial distribution. It implies that for any level of noise, the stiffness and the force of tension cannot be estimated. This observation is consistent with the expected results. There is indeed an infinity of solutions $(\alpha, f)$ that can verify $u(\alpha, f) = u_{\text{meas}}$, as long as the relation $f = \alpha u_{\text{meas}}$ is respected. It is therefore impossible to identify both parameters simultaneously.

<figure><img src="image_19.png" /><figcaption>Figure 15: Error distribution on the estimated stiffness (left) and on the estimated force (right) when both are estimated simultaneously</figcaption></figure>

The results presented so far do not take into account the mass of the spring. However, in order to validate our method, which consists in taking gravity into account to improve the estimation quality, we will include the mass of the spring from now on. Since the method is applied to the lung, we can draw a parallel between a lung and a spring.

Let us use the notation $+$ for all the expressions associated with the case of the lung in prone position, and $-$ for all the expressions associated with the case in supine position. The displacement of the spring associated to the prone and supine positions can therefore be written:

$$ u^+ = \frac{f + mg}{k} \qquad (46) $$

and

$$ u^- = \frac{f - mg}{k}. \qquad (47) $$

Let us try to identify simultaneously the stiffness and the force of tension of the spring. The cost function of the prone and supine positions can be written:

$$ J^+ = \frac{\|u^+(\theta) - u_{\text{meas}}^2\|^2}{\|n\|^2} \qquad (48) $$

and

$$ J^- = \frac{\|u^-(\theta) - u_{\text{meas}}^2\|^2}{\|n\|^2} \qquad (49) $$

Our new case, combining images in prone and supine positions, uses the following cost function, combining the two cost functions previously presented:

$$ J = \frac{\|u^+(\theta) - u_{\text{meas}}^+\|^2}{\|n\|^2} + \frac{\|u^-(\theta) - u_{\text{meas}}^2\|^2}{\|n\|^2} \qquad (50) $$

The distribution of the error on the estimated stiffness parameter and on the estimated force of tension, when both parameters are estimated simultaneously, are similar to the distributions
presented in Figure 15 for the prone and supine positions. As presented previously, it is indeed impossible to identify both parameters simultaneously, since there is an infinity of $(\alpha, p)$ pairs allowing us to retrieve the measured displacement field. However, combining images in prone and supine positions allows us to identify simultaneously both parameters.

The distribution of the error on the estimated stiffness parameter and on the estimated force of tension, when both parameters are estimated simultaneously, for the case combining prone and supine images, is presented in Figure 16. For any level of noise, our method provides the estimation error for the stiffness and the force of tensions. For large SNRs, the optimization always converges towards the exact parameter value; whereas for small SNRs, the obtained parameters are far away from the ground truth. Both parameters are therefore identifiable for reasonable SNRs when combining images in prone and supine positions, which is not possible when considering images only in prone or only in supine position. Implementing gravity in our model is therefore very interesting in our case, since it could allow us to identify additional parameters.

<figure><img src="image_20.png" /><figcaption>Figure 16: Error distribution on the estimated stiffness (left) and on the estimated force (right) when both are estimated simultaneously, when combining images in prone and supine positions.</figcaption></figure>
## Appendix C Recall of the weak formulation of the problem

This appendix details the weak formulation of the problem solved in this article. Details can be

found in [Peyraut et al. 2024b].

Based on the developments presented in the Methods, we can write the weak formulation of the

unloading problem as

$$ (\underline{u}, \phi_{s0}) \mid \left\{ \begin{array}{l} \displaystyle \int_{\omega} \underline{\sigma}_s : \underline{\epsilon}(\underline{u}^*) d\omega - \int_{\omega} p_f \operatorname{tr}(\underline{\epsilon}(\underline{u}^*)) d\omega = \int_{\omega} \rho \phi_s \underline{g} \cdot \underline{u}^* d\omega + \int_{\gamma} p_{pl} \underline{n} \cdot \underline{u}^* d\gamma \quad \forall \underline{u}^* \\ \\ \displaystyle \int_{\omega} (p_f + \frac{\partial \Psi_s}{\partial \Phi_s}) \phi_{s0}^* d\omega = 0 \quad \forall \phi_{s0}^* \end{array} \right. , \quad (51) $$

with $p_{pl}$ defined through the following system

$$ (p_{pl}, \underline{\lambda}, \underline{\mu}, \delta) \mid \left\{ \begin{array}{l} \displaystyle \int_{\gamma} (\underline{\nabla}_s p_{pl} - \underline{\nabla}_s \tilde{p}_{pl}) \cdot \underline{\nabla}_s p^* d\gamma - \underline{\lambda} \cdot \int_{\gamma} p^* \underline{n} d\gamma \\ \qquad + \underline{\mu} \cdot \int_{\gamma} (\underline{x} - \tilde{\underline{x}}) \times p^* \underline{n} d\gamma + \delta \int_{\gamma} p^* d\gamma = 0 \quad \forall p^* \\ \\ \displaystyle \int_{\omega} \rho \phi_s \underline{g} d\omega - \int_{\gamma} p_{pl} \underline{n} d\gamma = 0 \\ \\ \displaystyle \int_{\gamma} (\underline{x} - \tilde{\underline{x}}) \times p_{pl} \underline{n} d\gamma = 0 \\ \\ \displaystyle \int_{\gamma} (p_{pl} - \tilde{p}_{pl}) d\gamma = 0 \end{array} \right. . \quad (52) $$

Similarly, the weak formulation of the loading problem is expressed as

$$ (\underline{U}, \underline{\Phi}_s) \mid \left\{ \begin{array}{l} \displaystyle \int_{\Omega} \underline{\underline{\Sigma}}_s : d_{\underline{U}} \underline{\underline{E}} \cdot \underline{U}^* d\Omega - \int_{\Omega} p_f J \underline{\underline{C}}^{-1} : d_{\underline{U}} \underline{\underline{E}} \cdot \underline{U}^* d\Omega = \int_{\Omega} \rho \underline{\Phi}_{s0} \underline{g} \cdot \underline{U}^* d\Omega \\ \\ \displaystyle \qquad + \int_{\Gamma} p_{pl} J \underline{\underline{\Gamma}}^{-T} \cdot \underline{\underline{N}} \cdot \underline{U}^* d\Gamma \quad \forall \underline{U}^* \\ \\ \displaystyle \int_{\Omega} (p_f + \frac{\partial \Psi_s}{\partial \Phi_s}) \underline{\Phi}_s^* d\Omega = 0 \quad \forall \underline{\Phi}_s^* \end{array} \right. , \quad (53) $$

with $p_{pl}$ defined as in Equation (52).