# A quasi-static poromechanical model of the lungs

Cécile Patte, Martin Genet, Dominique Chapelle

► To cite this version:

Cécile Patte, Martin Genet, Dominique Chapelle. A quasi-static poromechanical model of the lungs. Biomechanics and Modeling in Mechanobiology, 2022, 21 (2), pp.527-551. 10.1007/s10237-021-01547-0. hal-03474200v2

HAL Id: hal-03474200

https://inria.hal.science/hal-03474200v2

Submitted on 17 Jan 2022

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
# A quasi-static poromechanical model of the lungs

Cécile Patte¹,², Martin Genet²,¹, Dominique Chapelle¹,²ᵃ

¹ Inria, {cecile.patte, dominique.chapelle}@inria.fr

² LMS, École Polytechnique, Institut Polytechnique de Paris, CNRS, martin.genet@polytechnique.edu

ᵃCorresponding author

**Abstract** — The lung vital function of providing oxygen to the body heavily relies on its mechanical behavior, and the interaction with its complex environment. In particular, the large compliance and the porosity of the pulmonary tissue are critical for lung inflation and air inhalation, and the diaphragm, the pleura, the rib cage and intercostal muscles all play a role in delivering and controlling the breathing driving forces. In this paper, we introduce a novel poromechanical model of the lungs. The constitutive law is derived within a general poromechanics theory via the formulation of lung-specific assumptions, leading to a hyperelastic potential reproducing the volume response of the pulmonary mixture to a change of pressure. Moreover, physiological boundary conditions are formulated to account for the interaction of the lungs with their surroundings, including a following pressure and bilateral frictionless contact. A strategy is established to estimate the unloaded configuration from a given loaded state, with a particular focus on ensuring a positive porosity. Finally, we illustrate through several realistic examples the relevance of our model and its potential clinical applications.

**Keywords** — Pulmonary Mechanics, Modeling, Poromechanics, Inverse Poromechanics, Finite Element Method.

## 1 Introduction

**Context and motivations** Pulmonary diseases are among the main causes of death in the world [OMS, 2016] and represent an important health burden. Among them, interstitial lung diseases (ILD) affect the pulmonary tissue and the organ structure, and consequently alter the pulmonary function. Indeed, the main lung function, consisting in providing oxygen to the whole body, relies on the alveoli structure of the lungs. The large compliance of the lungs throughout the breathing cycle also plays a role, letting air enter the bronchial tree all the way to the alveoli. The lung structure, mechanics and function are thus tightly linked. Moreover, mechanics is also assumed to be a key factor in the development of some diseases like pulmonary fibrosis, in which a vicious circle, where fibrosis leads to large stresses which in turn induce fibrosis spread, may occur [Hinz and Suki, 2016].

Understanding lung physiology allows to propose mathematical model equations for the phenomena at hand, and develop associated computational approaches. Computational medicine is a rapidly growing field benefitting from the progress in mathematical modeling and computational tools. Computational modeling is well developed in some medical disciplines like cardiology where actual clinical applications have already been addressed [Winslow et al., 2012; Lee et al., 2014; Chabiniok et al., 2016], whereas numerical pulmonology is still an emergent area [Tawhai et al., 2009; Kaul, 2019; Patte et al., 2020]. Early works in lung computational modeling have been mainly
motivated by image registration for radiotherapy purposes [Eom et al., 2010], and by airflow study in asthma or chronic obstructive pulmonary disease (COPD) [Bordas et al., 2015]. High hopes are placed in these new computational tools designed for clinical applications. They pave the way for improvements of diagnosis and prognosis, and can be used to test various therapeutic strategies in silico, thereby enabling the reduction of risks for patients and healthcare costs, as well as improved personalized care for each patient.

**Lung mechanical modeling** Lung modeling, and especially pulmonary mechanical modeling, is a great challenge given the high level of complexity of the organ and its environment [Clark et al., 2017]. A major difficulty lies in the many scales involved in pulmonary physiology and pathology, from the organ to the alveoli, through the lobe, segment, lobule and acinus scales [Burrowes et al., 2019].

A number of lung mechanical models have been proposed directly at the organ scale, with varying complexity depending on their objectives. Radiotherapy applications mainly require the modeling of displacements, that can be acquired in part with image registration [Sundaram and Gee, 2005; Hurtado et al., 2017]. More advanced mechanical models integrate more complex modeling ingredients for boundary conditions and material constitutive behavior. Some models apply a given displacement field on the lung boundary – or an orthogonal projection thereof – obtained by image registration [Brock et al., 2005; Al-Mayah et al., 2011]. Tawhai et al. [2009] define the so-called pleural cavity surface, along which the lung surface undergoes frictionless sliding. Other models apply a negative pressure representing the pleural pressure on the lung surface [Baudet et al., 2003; Zhang et al., 2004; Werner et al., 2009; Han et al., 2017], while still prescribing the final shape by means of a contact surface.

As for material constitutive behavior, lung parenchyma is often modeled with a linear isotropic stress-strain law, either within a totally linear formulation [Baudet et al., 2003; Zhang et al., 2004; Brock et al., 2005; Werner et al., 2009; Fuerst et al., 2015], or with a neo-Hookean hyperelastic law [Berger et al., 2016; Han et al., 2017]. In order to take into account the non-linearity of the lung behavior, West and Matthews [1972] use a strain-dependent value of Young's modulus. Zeng et al. [1987] propose a more complex hyperelastic potential based on experimental mechanical tests to capture the non-linear behavior, and this hyperelastic potential was subsequently used by [Eom et al., 2010; Al-Mayah et al., 2011]. However, as it has been formulated based on biaxial tension tests, this constitutive law cannot adequately reproduce the pressure-volume behavior measured in inflation tests such as those of Richardson et al. [2019]. By contrast, Tawhai et al. [2009] use a hyperelastic potential able to reproduce the pressure-volume relationship.

Another key modeling aspect concerns the treatment of the unloaded configuration in the mechanical formulation, as the lung tissue is a pre-stressed system that operates under varying tensile stresses – balancing the negative pleural pressure – over the whole breathing cycle. Most studies using a non-linear behavior either neglect the fact that the initial configuration extracted from images is not stress-free [Eom et al., 2010; Al-Mayah et al., 2011], or use a direct scaling to get the supposedly unloaded configuration from the initial one [Tawhai et al., 2009; Berger et al., 2016].

Other modeling approaches have been considering the lung microstructure, in order to investigate structure-properties relationships. Models of the alveoli microstructure have been developed using given geometrical shapes [Dale et al., 1980; Budiansky and Kimmel, 1987; Denny and Schroter, 2006], or image segmentation [Sarabia-Vallejos et al., 2019; Álvarez-Barrientos et al., 2021]. The link between organ scale and alveolar scale can then be performed by homogenization techniques [Wiechert and Wall, 2010; Cazeaux and Grandmont, 2015], but this is a computationally costly
approach especially when non-linear behavior is considered.

An alternative approach to take into account the coupling of air flows and parenchyma without explicitly modeling the microstructure is provided by poromechanics. General poromechanics theories [Biot and Temple, 1972; Coussy, 2004; Dormieux et al., 2006] have been widely used in geophysics, and have then been adapted for the needs of biomechanics by Chapelle and Moireau [2014] with a theory compatible with large strains and rapid flows. Such a poromechanical framework has been successfully applied in biomechanics for cardiac perfusion [Chapelle et al., 2010; Reeve et al., 2014; Burtschell, 2016] or bones [Scheiner et al., 2013]. For lungs, beyond the preliminary work of Kowalczyk [1993], poromechanics approaches are quite scarce, and mainly focused on air flows [Seyfi et al., 2016; Berger et al., 2016].

**Our work** In this work, we propose to model lung mechanical behavior using poromechanics. We use the general poromechanics theory of Chapelle and Moireau [2014], which is compatible with large deformations, and we adapt it to the pulmonary case with the incorporation of lung-specific assumptions. Our model includes physiological boundary conditions and a hyperelastic potential adequately reproducing the volume response of the pulmonary mixture to a change of pressure. Our ultimate objective is to perform patient-specific modeling, and therefore our complete modeling framework is designed to be amenable to that purpose. A strategy is presented to estimate the unloaded configuration with a particular focus on ensuring a positive porosity, since this is an actual issue when large displacements are considered. Finally, we provide several illustrations of the ability of our model to represent various types of physiological or pathological conditions, paving the way for future clinical applications.

# 2 Lung model

## 2.1 Poromechanical framework

With the motivations given in the introduction, we propose in this paper a model of lung biomechanics at the organ space scale and the breathing time scale. The model is written in the general framework of Biot's macroscopic theory of poromechanics [Biot, 1941; Biot and Temple, 1972; Coussy, 2004; Chapelle and Moireau, 2014]. Thus, the lung constituents (i.e. the parenchymal tissue, air and blood) are not represented explicitly and separately, but via the so-called mixture theory with two phases in which “solid” (the associated quantities of which being denoted with the subscript “s”) and “fluid” (subscript “f”) are distinct distributed phases present everywhere in varying fractions and interacting with each other. We choose to consider the air in the airways and alveoli as the so-called fluid phase, and the rest, i.e. tissue and blood, as the solid phase.

### 2.1.1 Basic assumptions for the pulmonary setting

Beside standard considerations in continuum poromechanics [Coussy, 2004], the proposed model is based on the following specific assumptions:

* The mixture temperature is constant and uniform, i.e. body temperature. This is a natural hypothesis for biological tissues.

* The fluid is incompressible, as the Reynolds number of the airflows inside the lungs is low [Baffico et al., 2010].
* The breathing stages considered are associated with static equilibrium, which is consistent with clinical data that can be obtained with patients, see [Gibson and Pride, 1976]. Inertia and viscous terms are then neglected, both for the fluid and the solid phases.

* The fluid pressure is homogeneous throughout the lungs, as a consequence of the previous assumption. Its actual value depends on the breathing configuration considered: in the case of free-breathing the fluid pressure is equal to the atmospheric pressure, whereas in the case of ventilated breathing it is prescribed by the ventilator.

The first two assumptions have already been used in [Chapelle and Moireau, 2014] and are also valid in the case of cardiac perfusion modeling for example [Burtschell, 2016]. The last two are specific to our setting of pulmonary modeling. Their validity depends on the actual breathing conditions, and they are all the more valid as the breathing is slowly-varying. As a consequence, our model applies to breath-holding as well as slowly-varying phases such as end-exhalation and end-inhalation in normal breathing, but not to fast breathing.

### 2.1.2 Continuum poromechanics

**Kinematics** We consider a mixture deforming under some loading, and occupying the domain $\omega$ with boundary $\gamma$. Following the total Lagrangian approach, we denote by $\Omega_0$ a supposedly known reference configuration (which is not necessarily unloaded or stress-free, or even equal to $\omega_{t=0}$) with boundary $\Gamma_0$, on which we will formulate the problem, so that all quantities must be defined or transported onto $\Omega_0$. Quantities defined over the reference configuration will be written in uppercase, while quantities defined over the deformed configuration will be written in lowercase. Furthermore, to distinguish quantities characterizing the reference and the current configurations, we will write the former with a subscript “0”. The deformation, i.e. the mapping from $\Omega_0$ to $\omega$, is denoted by $\chi$, and the displacement is given by

$$ \underline{U}(\underline{X}) := \underline{\chi}(\underline{X}) - \underline{X} = \underline{x}(\underline{X}) - \underline{X}. \qquad (1) $$

The deformation gradient is then

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

so that the relative volume change of the mixture is

$$ J := \det \underline{F}. \qquad (3) $$

We also introduce the right Cauchy-Green deformation tensor

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

as well as the Green-Lagrange strain tensor

$$ \underline{E} := \frac{1}{2} (\underline{C} - \underline{1}). \qquad (5) $$
<figure><img src="image_2.png" /><figcaption>Figure 1: Schematic of the reference and deformed configurations with associated local quantities, volumes and porosity. $d\Omega_0$, $d\omega$, $d\Omega_{f0}$, $d\omega_f$, $d\Omega_{s0}$ $d\omega_s$ are the infinitesimal volumes of the mixture, the fluid and solid phases, in the reference and deformed configurations, respectively. $\Phi_{f0}$ and $\phi_f$ are the porosities in the reference and deformed configurations, respectively. $\bar{\rho}_{f+}$ is the added fluid mass density defined in the reference configuration. $J$, $\Phi_s$, $\Phi_f$ are the local volume ratios.</figcaption></figure>
**Specific poromechanics variables** To formally define poromechanics variables, we consider the infinitesimal domains of matter surrounding X and x, denoted by $d\Omega_0$ and $d\omega$, respectively. Both can be decomposed into solid and fluid parts, i.e. $d\Omega_0 = d\Omega_{f0} \cup d\Omega_{s0}$ and $d\omega = d\omega_f \cup d\omega_s$. It is important to note that $d\omega$ is the image of $d\Omega_0$ through the mapping $\chi$, so that it contains the same solid matter; the fluid, however, can flow in and out of $d\omega$ and is not necessarily conserved locally.

The key quantity of poromechanics is the porosity, defined as the volume fraction of fluid in the mixture. We will denote

$$ \Phi_{f0} := \frac{|d\Omega_{f0}|}{|d\Omega_0|} \quad \text{and} \quad \phi_f := \frac{|d\omega_f|}{|d\omega|} \qquad (6) $$

the porosity in the reference and the current configurations, respectively.

We also introduce

$$ \Phi_s := J \left(1 - \phi_f \circ \chi\right) = \frac{|d\omega_s|}{|d\Omega_0|} \qquad (7) $$

the current volume fraction of solid pulled back on the reference configuration, and we denote by

$$ \Phi_f := J \phi_f \circ \chi = \frac{|d\omega_f|}{|d\Omega_0|} \qquad (8) $$

the current volume fraction of fluid pulled back on the reference configuration, such that $\Phi_f = J - \Phi_s$.

We now introduce mass densities, distinguishing between constituent (fluid or solid) actual densities, which are defined as mass per unit constituent volume, from apparent densities, which are defined as mass per unit mixture volume and written with a bar. Thus, the reference fluid and solid mass densities are denoted by $\rho_{f0}$ and $\rho_{s0}$, respectively. The mixture reference mass density is decomposed into a fluid and a solid part

$$ \bar{\rho}_0 = \bar{\rho}_{f0} + \bar{\rho}_{s0}, \qquad (9) $$

with $\bar{\rho}_{f0} := \rho_{f0}\Phi_{f0}$ and $\bar{\rho}_{s0} := \rho_{s0}(1 - \Phi_{f0})$ the fluid and solid mass per unit volume of the mixture in the reference configuration, i.e. the reference apparent densities.

In the current configuration we denote the fluid and solid mass densities by $\rho_f$ and $\rho_s$, respectively, and we also introduce the associated quantities pulled back to the reference configuration, i.e. $\bar{\rho}_f := \rho_f J \phi_f \circ \chi$ and $\bar{\rho}_s := \rho_s J (1 - \phi_f \circ \chi)$, leading to the total pulled back apparent density

$$ \bar{\rho} = \bar{\rho}_f + \bar{\rho}_s. \qquad (10) $$

Note that $\bar{\rho}_f$, $\bar{\rho}_s$ and $\bar{\rho}$ represent current mass per unit reference mixture volume, and that solid mass conservation entails

$$ \bar{\rho}_s = \bar{\rho}_{s0}. \qquad (11) $$

We can now define the “added” fluid mass per unit reference mixture volume as

$$ \bar{\rho}_{f+} := \bar{\rho}_f - \bar{\rho}_{f0} = \bar{\rho} - \bar{\rho}_0. \qquad (12) $$

Using the fluid incompressibility assumption whereby $\rho_f = \rho_{f0}$, we get

$$ \bar{\rho}_{f+} = \rho_{f0} (\Phi_f - \Phi_{f0}). \qquad (13) $$

This added fluid mass quantity can be positive (inflation) or negative (deflation). However, one cannot remove more fluid than there initially was, so that $\bar{\rho}_{f+} \ge -\bar{\rho}_{f0}$.
**Force balance** At any time the mixture is in equilibrium with external loadings, such as boundary tractions or gravity. The stress tensor is denoted by $\underline{\underline{\sigma}}$ in $\omega$ (Cauchy stress tensor) and $\underline{\underline{\Sigma}} := J\underline{\underline{F}}^{-1} \cdot \underline{\underline{\sigma}} \circ \underline{\underline{\chi}} \cdot \underline{\underline{F}}^{-T}$ in $\Omega_0$ (second Piola-Kirchhoff stress tensor). Mechanical equilibrium of the mixture can be expressed in weak form through the static principle of virtual work written in the reference configuration as

$$ \forall \underline{U}^*, W_{\text{int}}(\underline{U}, \underline{U}^*) = W_{\text{ext}}(\underline{U}, \underline{U}^*), \quad (14) $$

where $W_{\text{int}}$ is the virtual work – associated with the virtual displacement field $\underline{U}^*$ – of internal forces of the mixture and $W_{\text{ext}}$ is the virtual work of external forces applied to the mixture. The virtual work of internal forces can be expressed as

$$ W_{\text{int}}(\underline{U}, \underline{U}^*) = \int_{\Omega_0} \underline{\underline{\Sigma}} : d\underline{U} \underline{\underline{E}} \cdot \underline{U}^* d\Omega_0, \quad (15) $$

where $d\underline{U} \underline{\underline{E}} \cdot \underline{U}^* = (\underline{\underline{F}}^T \cdot \underline{\underline{\nabla}} \underline{U}^*)_{\text{sym}}$ is the differential of the Green-Lagrange strain tensor. The virtual work of internal forces $W_{\text{int}}$ will be detailed in Section 2.1.3 through the choice of a free energy function, while the virtual work of external forces $W_{\text{ext}}$ will be detailed in Section 2.3 through the definition of loading conditions. General poromechanics formulations usually require other balance equations for the fluid flow and fluid mass conservation [Chapelle and Moireau, 2014]. In our specific setting, the static assumption makes the fluid flow equation trivial, and the fluid mass conservation is not designed to be exploited in a static configuration, but it is made superfluous by the assumption of prescribed fluid pressure that allows to close the system, as will be seen in the next section.

### 2.1.3 Constitutive behavior

**Mixture free energy** The mixture constitutive behavior is characterized through its Helmholtz free energy density potential denoted by $\bar{\psi}$, which is a function of the Green-Lagrange strain tensor $\underline{\underline{E}}$ and added fluid mass $\bar{\rho}_{f+}$ [Coussy, 2004; Dormieux et al., 2006; Chapelle and Moireau, 2014]. The fundamental principle of mixture theories states that it is additively decomposed into solid and fluid parts as

$$ \bar{\psi} = \bar{\psi}_s + \bar{\psi}_f, \quad (16) $$

where $\bar{\psi}_s$ and $\bar{\psi}_f$ are the solid and fluid free energies per unit mixture volume in the reference configuration.

It is convenient to express the solid free energy density $\bar{\psi}_s$ as a function of $\underline{\underline{E}}$ and of the pulled back solid volume fraction $\Phi_s$. Indeed, we will see that the energy conjugate of $\Phi_s$ is directly involved in the fundamental poromechanical equilibrium principle between fluid and solid hydrostatic stresses, making it a very natural variable, and this choice is possible since $\Phi_s$ can be expressed from the state variables $\underline{\underline{E}}$ and $\bar{\rho}_{f+}$, as seen from Eq. (13) and $\Phi_s = J - \Phi_f$.

Regarding the fluid, since it is assumed to be incompressible and in isothermal conditions, its free energy density can be expressed as a function of $\Phi_f$ only, i.e.

$$ \bar{\psi}_f(\Phi_f) = -p_{f0} \Phi_f, \quad (17) $$

where $p_{f0}$ is a reference pressure, see [Chapelle and Moireau, 2014] for more detail. Note that $\Phi_f$ depends only on the state variable $\bar{\rho}_{f+}$, recall Eq. (13).
Finally, we consider a mixture free energy density in the form:

$$ \bar{\psi}(\underline{E}, \bar{\rho}_{f+}) = \bar{\psi}_s(\underline{E}, \Phi_s) + \bar{\psi}_f(\Phi_f), \quad (18) $$

with the identities

$$ \left\{ \begin{aligned} \Phi_s &= J(\underline{E}) - \Phi_f \\ \Phi_f &= \Phi_{f0} + \frac{\bar{\rho}_{f+}}{\rho_{f0}} \end{aligned} \right. \quad (19) $$

Since all dissipative processes have been neglected, the second principle of thermodynamics entails that the second Piola-Kirchhoff stress tensor simply derives from the free energy potential. We then have

$$ \begin{aligned} \underline{\Sigma} &= \left. \frac{\partial \bar{\psi}(\underline{E}, \bar{\rho}_{f+})}{\partial \underline{E}} \right|_{\bar{\rho}_{f+}} = \left. \frac{\partial \bar{\psi}_s(\underline{E}, \Phi_s)}{\partial \underline{E}} \right|_{\bar{\rho}_{f+}} + \left. \frac{\partial \bar{\psi}_f(\Phi_f)}{\partial \underline{E}} \right|_{\bar{\rho}_{f+}} \\ &= \left. \frac{\partial \bar{\psi}_s(\underline{E}, \Phi_s)}{\partial \underline{E}} \right|_{\Phi_s} + \left. \frac{\partial \bar{\psi}_s(\underline{E}, \Phi_s)}{\partial \Phi_s} \right|_{\underline{E}} \left. \frac{\partial \Phi_s}{\partial \underline{E}} \right|_{\bar{\rho}_{f+}} \\ &= \frac{\partial \bar{\psi}_s(\underline{E}, \Phi_s)}{\partial \underline{E}} - p_s J \underline{C}^{-1}, \end{aligned} \quad (20) $$

where $p_s := -\frac{\partial \bar{\psi}_s(\underline{E}, \Phi_s)}{\partial \Phi_s}$ is the solid hydrostatic stress induced by volume fraction changes. Internal equilibrium of the mixture requires that [Coussy, 2004; Chapelle and Moireau, 2014]

$$ p_s = p_f, \quad (21) $$

leading to

$$ \underline{\Sigma} = \frac{\partial \bar{\psi}_s(\underline{E}, \Phi_s)}{\partial \underline{E}} - p_f J \underline{C}^{-1}, \quad (22) $$

together with

$$ p_f = - \frac{\partial \bar{\psi}_s(\underline{E}, \Phi_s)}{\partial \Phi_s}. \quad (23) $$

We will now define a specific solid free energy density for the parenchyma.

**Solid free energy** Following [Chapelle and Moireau, 2014], we choose an additive form of the solid free energy density potential split into skeleton and bulk parts, i.e.

$$ \bar{\psi}_s(\underline{E}, \Phi_s) = \bar{W}_{\text{skel}}(\underline{E}) + \bar{W}_{\text{bulk}}(\Phi_s), \quad (24) $$

where $\bar{W}_{\text{skel}}$ represents the mechanical behavior of the porous solid, or skeleton, as a structure subjected to deformation, and $\bar{W}_{\text{bulk}}$ its mechanical behavior with regard to solid volume fraction changes. Such a decomposition leads to the following expressions for the stress and internal hydrostatic equilibrium

$$ \underline{\Sigma} = \frac{\partial \bar{W}_{\text{skel}}(\underline{E})}{\partial \underline{E}} - p_f J \underline{C}^{-1} \quad \text{and} \quad p_f = - \frac{\partial \bar{W}_{\text{bulk}}(\Phi_s)}{\partial \Phi_s}. \quad (25) $$
It is possible to define these potentials directly, assuming a given solid volume fraction, in which case we will call them *effective potentials*, representing an energy density per unit volume of the mixture in the reference configuration. Alternatively, we can scale these energy densities by the solid volume fraction, i.e. introducing $\tilde{\psi}_s$, $\tilde{W}_{skel}$ and $\tilde{W}_{bulk}$ such that

$$ 
\left\{ 
\begin{aligned}
\tilde{\psi}_s &= (1 - \Phi_{f0}) \tilde{\psi}_s \\
\tilde{W}_{skel} &= (1 - \Phi_{f0}) \tilde{W}_{skel} \\
\tilde{W}_{bulk} &= (1 - \Phi_{f0}) \tilde{W}_{bulk}
\end{aligned}
\right. 
\qquad (26) $$

These rescaled quantities $\tilde{\psi}_s$, $\tilde{W}_{skel}$ and $\tilde{W}_{bulk}$ are now energy densities per unit solid volume in the reference configuration, and we will call them *rescaled potentials*. Note that we do not expect these rescaled potentials to be independent of the porosity, since structure effects – associated with the microstructure – intervene in a manner that is not directly proportional to mass. Nevertheless, mass effects should be dominant in the bulk energy, so we expect $\tilde{W}_{bulk}$ to have a weak dependence on porosity.

The pulmonary parenchyma is assumed to behave as an isotropic, non-linear, hyperelastic material. Even though it is assumed to be quasi-incompressible as a material, its response as a structure is not at all incompressible, since fluid can flow in and out of each local volume. The modeling of the quasi-incompressibility of the solid material will be dealt with in the definition of the bulk potential below, whereas we use the following compressible potential for the skeleton free energy

$$ \tilde{W}_{skel}(\underline{E}) = \tilde{\beta}_1 (I_1 - 3 - 2 \ln J) + \tilde{\beta}_2 (I_2 - 3 - 4 \ln J) + \tilde{\alpha} (e^{\delta(J^2 - 1 - 2 \ln J)} - 1), \quad (27) $$

where $I_1 := \text{tr } \underline{C}$ and $I_2 := \frac{1}{2}(\text{tr}(\underline{C})^2 - \text{tr}(\underline{C}^2))$ are the first two invariants of $\underline{C}$, and $\tilde{\beta}_1, \tilde{\beta}_2, \tilde{\alpha}$ and $\delta$ are four material parameters. This potential contains standard neo-Hookean and Mooney-Rivlin terms [Rivlin, 1948; Ogden, 1972], and an exponential of the Ciarlet-Geymonat bulk term [Ciarlet and Geymonat, 1982]. The use of an exponential term, initiated by Demiray [1972] and popularized by Fung [1981], is common in soft tissues modeling [Weiss et al., 1996; Caruel et al., 2014; Genet et al., 2015a]. Here it was motivated by the shape of the clinical data representing the lung volume with respect to the pleural pressure [Gibson and Pride, 1976].

The effective hyperelastic potential $\tilde{W}_{skel}$ then depends on the four material constants $\tilde{\alpha}, \tilde{\beta}_1, \tilde{\beta}_2, \delta$, with the first three defined as

$$ 
\left\{ 
\begin{aligned}
\tilde{\alpha} &= (1 - \Phi_{f0}) \tilde{\alpha} \\
\tilde{\beta}_1 &= (1 - \Phi_{f0}) \tilde{\beta}_1 \\
\tilde{\beta}_2 &= (1 - \Phi_{f0}) \tilde{\beta}_2
\end{aligned}
\right. 
\qquad (28) $$

the typical values of which can be determined by calibration or identification based on clinical or experimental data. The volume macroscopic behavior of the pulmonary mixture is usually characterized by P-V curves, with P the pleural pressure and V the lung volume, that are used by clinicians to assess lung compliance [Gibson and Pride, 1976]. These data then allow to characterize the mixture volume changes with the pressure and to determine the dependency of the skeleton deformation energy $\tilde{W}_{skel}$ with $J$. We determined a set of the four effective parameters of the skeleton hyperelastic potential compatible with the data from [Gibson and Pride, 1976] over a range of physiological volumes, see Figure 2. As regards the well-posedness of the corresponding identification problem, it would clearly deserve further investigation. In our case, however, the
role of the neo-Hookean and Mooney-Rivlin terms is essentially to regularize the behavior near the reference configuration since the exponential term itself gives a very soft behavior for small deformations, and therefore we did not attempt to fine-tune the parameters $\beta_1$ and $\beta_2$. Typical values for a healthy case were obtained as $\bar{\alpha} = 0.08$ kPa, $\bar{\beta}_1 = 0.1$ kPa, $\bar{\beta}_2 = 0.2$ kPa, $\delta = 0.5$.

For the bulk term $\overline{W}_{\text{bulk}}$, we use the following expression, as proposed in [Chapelle and Moireau, 2014]:

$$ \tilde{W}_{\text{bulk}}(\Phi_s) = \tilde{\kappa} \left( \frac{\Phi_s}{(1 - \Phi_{f0})} - 1 - \ln \left( \frac{\Phi_s}{(1 - \Phi_{f0})} \right) \right), \quad (29) $$

where $\tilde{\kappa}$ is the bulk modulus of the solid material, here assumed to be large to represent quasi-incompressibility of the solid. This potential is an extension of the Ciarlet-Geymonat bulk potential [Ciarlet and Geymonat, 1982] to poromechanics; however, other penalization terms with similar properties could have been used. Note also that – as is usual in solid mechanics – we prefer to consider quasi-incompressibility rather than exact incompressibility because this is more realistic and does not induce any serious difficulty [Le Tallec, 1994].

### 2.1.4 Methods to ensure positive porosity

The above poromechanical model, with the choices made on the energies in Section 2.1.3, does not ensure by itself a positive porosity [Chapelle et al., 2010; Chapelle and Moireau, 2014], which may be an issue with large strains during mixture compression or, in the pulmonary case, drastic deflation such as that associated with pneumothorax. We now discuss two methods to address this issue.

**Use of an additional energy** Chapelle and Moireau [2014] proposed adding a penalization term to the solid energy $\bar{\psi}_s$ in order to keep the porosity positive. The additional energy tends to infinity when the porosity $\phi_f$ tends to 0, which means that it would require an infinite amount of energy – or equivalently an infinitely negative fluid pressure – to close a pore completely, as is natural from a physical standpoint.

Denoting this additional term by $\overline{W}_{\text{por}}(\Phi_f; \Phi_{f0})$ – where $\Phi_{f0}$ appears as a parameter – the following specifications should be enforced:

* $\overline{W}_{\text{por}}(\Phi_f; \Phi_{f0})$ tends to infinity when $\Phi_f$ tends to zero (for any $\Phi_{f0}$), to serve the above-described purpose;

* $\overline{W}_{\text{por}}(\Phi_f; \Phi_{f0}) = 0$ for $\Phi_f/\Phi_{f0}$ in a wide range $[r_{\text{inf}}, r_{\text{sup}}]$ with $r_{\text{inf}} \ll 1$ and $r_{\text{sup}} \gg 1$, so that the corresponding contribution vanishes for a standard regime of behavior;

* $\overline{W}_{\text{por}}(\Phi_f; \Phi_{f0})$ tends to infinity when $\Phi_{f0}$ tends to zero for $\Phi_f$ fixed, which will be useful for the inverse problem of seeking the unloaded configuration with unknown porosity, see Section 2.5.4;

* $\overline{W}_{\text{por}}(\Phi_f; \Phi_{f0})$ of regularity $C^1$, and preferably convex for minimization purposes.

Note that the above specifications can be easily fulfilled by constructing $\overline{W}_{\text{por}}(\Phi_f; \Phi_{f0})$ based on a function of $r_\Phi = \Phi_f/\Phi_{f0}$, for instance

$$ \overline{W}_{\text{por}}(\Phi_f; \Phi_{f0}) := \begin{cases} \eta \left( \frac{r_{\text{inf}}}{r_\Phi} - 1 \right)^{n+1} & \text{if } r_\Phi < r_{\text{inf}} \\ 0 & \text{if } r_{\text{inf}} \le r_\Phi < r_{\text{sup}} \\ \eta \left( \frac{r_\Phi}{r_{\text{sup}}} - 1 \right)^{n+1} & \text{if } r_{\text{sup}} \le r_\Phi \end{cases} \quad (30) $$
<figure><img src="image_3.png" /><figcaption>Figure 2: Mechanical response of the mixture resulting from our choice of $\tilde{W}_{\text{skel}}$ in Eq. (27) & $\tilde{W}_{\text{bulk}}$ in Eq. (29) for various values of bulk modulus $\tilde{\kappa}$ (while other parameters are fixed at $\Phi_{f0} = 0.5$, $\tilde{\alpha} = 0.16$ kPa, $\tilde{\beta}_1 = 0.2$ kPa, $\tilde{\beta}_2 = 0.4$ kPa & $\delta = 0.5$) and applied pressure – positive when internal or negative when external, here considered in absolute value in the plots. The solid black curve corresponds to an external negative pressure loading (i.e., by imposing $\Sigma > 0$ and $p_f = 0$, in which case the response only depends on $\tilde{W}_{\text{skel}}$, not $\tilde{W}_{\text{bulk}}$, hence not on $\tilde{\kappa}$), while the colored curves correspond to a positive internal pressure loading (i.e., by imposing $p_f > 0$ and $\Sigma = 0$, in which case the response depends on $\tilde{\kappa}$, so various values are considered). All curves are obtained by imposing a level of pressure and solving the nonlinear behavior law at one single material point – i.e., a Rivlin cube – through a Newton-Raphson procedure, in order to compute the associated mixture deformation & volumes fractions changes. (Left) Volumetric response of a mixture to a pressure loading, compared to experimental data [Gibson and Pride, 1976]. $J_{\text{norm}}$ is the mixture volume change $J$ normalized by its value at 0.5 kPa (which roughly corresponds to the end-exhalation state in the lung model), leading to $J_{\text{norm}} = J/J_{0.5} = 1$ for an applied pressure of 0.5 kPa. The colored curves are superimposed with the black curve because the mixture volume deformation is identical for internal vs. external applied pressure. (Center & Right) Evolution of the volume fraction of fluid (i.e., the porosity) and solid with the pressure loading. When the mixture is subjected to internal pressure, the solid phase is compressed according to its bulk modulus; the lower the bulk modulus, the more compressible the solid, the higher the volume fraction of fluid and the lower the volume fraction of solid. For the value of the bulk modulus used in the rest of the paper ($\tilde{\kappa} = 10^2$ kPa) the change in solid volume is very limited.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>(a) Influence of the solid bulk modulus $\tilde{\kappa}$ when using the pore energy $\tilde{W}_{por}$. As long as the porosity is away from zero, the solid bulk modulus has little impact, since the volume decrease is mainly due to air ejection and less to solid compression. Then, when the porosity comes close to zero, the solid starts to be compressed according to the solid bulk modulus.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b) Influence of $\tilde{\eta}$ when using $\tilde{W}_{por}$. When $\tilde{\eta}$ decreases, the change of behavior when the porosity comes close to zero becomes sharper.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(c) Influence of the solid bulk modulus $\tilde{\kappa}$ when using the contact-like formulation. The same behavior as in Figure 3a is observed but the curves are sharper using the contact-like formulation.</figcaption></figure>

Figure 3: Behavior of a mixture under compression using either the pore energy formulation (here with $n = 4$ and $r_{inf} = r_{sup} = 1$) or the contact-like formulation. The mixture is in free-breathing condition ($p_f = 0$), and such that $\Phi_{f0} = 0.5$. In each case, the evolution of the mixture volume change $J$ and of the porosity $\Phi_f$ with respect to the pressure loading are shown. The black solid line illustrates the weakness of the model regarding negative porosity when neither the pore energy nor the contact-like formulation is used. The black dotted line represents $\Phi_f = 0$.
with $n > 1$ to ensure continuous differentiability. The material response induced by this potential is illustrated in Figure 3.

We thus consider

$$ \bar{\psi}_s(E, \Phi_s) = \bar{W}_{\text{skel}}(E) + \bar{W}_{\text{bulk}}(\Phi_s) + \bar{W}_{\text{por}}(\Phi_f), \quad (31) $$

with $\Phi_f = J - \Phi_s$, and (23) now yields

$$ p_f = - \frac{\partial \bar{W}_{\text{bulk}}(\Phi_s)}{\partial \Phi_s} + \frac{\partial \bar{W}_{\text{por}}(\Phi_f)}{\partial \Phi_f}, \quad (32) $$

which provides a new relationship between $p_f$, $\Phi_f$ and $J$, by which we see that $\Phi_f$ cannot tend to zero unless $p_f$ tends to minus infinity. Note also that the additional energy term is only a function of $\bar{\rho}_{f+}$ (via $\Phi_f$) and not of $E$, and therefore this does not change the expression of the stress tensor with respect to the energy densities, i.e.

$$ \Sigma = \frac{\partial \bar{\psi}(E, \bar{\rho}_{f+})}{\partial E} \bigg|_{\bar{\rho}_{f+}} = \frac{\partial \bar{W}_{\text{skel}}(E)}{\partial E} + \frac{\partial \bar{W}_{\text{bulk}}(\Phi_s)}{\partial \Phi_s} J C^{-1}. \quad (33) $$

However, the relationship between $\Sigma$ and $p_f$ is no longer as in (25) since the expression of $p_f$ itself has been modified.

**Contact-like formulation** When scaling the above penalization energy with a multiplicative coefficient $\eta$, the asymptotic case ($\eta \to 0$) is a contact-like problem, which can be formulated with two quantities, i.e. the porosity $\phi_f$ and the solid self-contact pressure $p_c$, using the following conditions

$$ \left\{ \begin{array}{l} \phi_f \cdot p_c = 0 \\ \phi_f \ge 0 \\ p_c \ge 0 \end{array} \right. \quad (34) $$

This expresses that there are two possible regimes for the porous material at every point: either it contains some fluid and the porosity is non-zero while the self-contact pressure is zero, or the pores are collapsed and the self-contact pressure is strictly positive while the porosity is zero. We will use the classical Heaviside step function $\mathcal{H}(\phi_f)$ such that

$$ \mathcal{H}(\phi_f) := \begin{cases} 1 & \text{if } \phi_f > 0 \\ 0 & \text{otherwise} \end{cases} \quad (35) $$

When $\mathcal{H}(\phi_f) = 1$, some fluid is present and the fluid pressure $p_f$ is well-defined, whereas $\mathcal{H}(\phi_f) = 0$ means that there is no fluid left so that fluid pressure is not a relevant quantity anymore. In the absence of fluid, the contact pressure $p_c$ is then a substitute of the fluid pressure and acts against the solid compression. We can define the aggregated pressure quantity $p$ as

$$ p := p_f \mathcal{H}(\phi_f) + p_c (1 - \mathcal{H}(\phi_f)), \quad (36) $$

which is equal to $p_f$ or $p_c$ depending on the state of the mixture. By extension of the internal equilibrium principle of Eq. (21) to the case $\phi_f = 0$, we have

$$ p = - \frac{\partial \bar{\psi}_s}{\partial \Phi_s}. \quad (37) $$
Initially, the porosity $\phi_f$ is non zero and the contact pressure $p_c$ is zero. When the mixture is compressed, the porosity decreases as fluid is expelled. Then, the porosity reaches 0, which means that the entire initial fluid has been expelled. If the mixture keeps being compressed, the solid part is compressed according to the solid bulk modulus $\kappa$, and the resulting contact pressure $p_c$ increases. This behavior is illustrated in Figure 3.

**Remark** As pointed out in [Chapelle and Moireau, 2014], the other constraint $\phi_f < 1$ is naturally satisfied with the proposed expression of the bulk energy $\overline{W}_{\text{bulk}}$ in Eq. (29). Indeed, when the porosity $\phi_f$ tends to 1, $\Phi_s$ tends to zero, and the bulk energy $\overline{W}_{\text{bulk}}$ tends to infinity.

## 2.2 Geometric model

Our lung model requires the definition of several geometric entities as illustrated in Figure 4, all described by finite element meshes in practice. In addition to the lungs themselves, which are volume meshes, two surface meshes are used, i.e. the rib cage surface and the inner surface. The rib cage surface is in contact with the lung where the lung is surrounded by hard tissue (i.e. ribs and spine). It corresponds to the costal pleura. The inner surface is in contact with the lung where the lung is surrounded by soft tissues. It corresponds to the diaphragmatic pleura and mediastinal pleura, i.e. the soft surroundings of lungs. Both rib cage and inner surfaces are surface meshes. The interactions between lungs and these surfaces are discussed in Section 2.3, and a method to produce these surfaces is presented in Section 2.6.

Typically, realistic geometric models can be obtained from biomedical images through segmentation. However, it is important to notice that no image of the unloaded configuration can be acquired *in vivo* as this would correspond to a pneumothorax, in which air enters the pleural space and pleural pressure becomes equal to atmospheric pressure. Consequently, the unloaded configuration associated with image-based geometries must be computed, as described in Section 2.4.

## 2.3 Loading and boundary conditions

**Elements of pulmonary physiology** A complex environment is involved during breathing. The lungs are surrounded by the pleura, composed of two membranes and applying a negative pressure, called pleural pressure, on the lungs surface. The pleural pressure keeps the lungs inflated and maintains the contact with their close environment (thorax and diaphragm mainly). Two mechanisms of breathing should be distinguished, spontaneous breathing and mechanical ventilation, differing by the driving force for lung inflation. During spontaneous breathing, the diaphragm contracts and pulls on the external layer of the pleura, causing a decrease in the pleural pressure, which induces lungs inflation and air filling. Intercostal muscles also play a role in the thorax inflation during inhalation. By contrast, during mechanical ventilation the diaphragm and intercostal muscles do not contract. The driving force is then provided by the ventilator, which induces a positive air pressure inside the lungs to make them inflate. This pressure is also present near the lung surface, where it pushes against the rib cage and the diaphragm, which themselves tend to resist lung inflation, contrary to spontaneous breathing.

**Spontaneous breathing** As shown in Figure 4, several choices are made to model the complex loading and boundary conditions for spontaneous breathing.
(i) The air pressure is homogeneous and equal to the atmospheric pressure according to the assumptions of Section 2.1.1, i.e. $p_f = 0$.

(ii) The pleural pressure $p_{pl}$ is modeled by a negative pressure $p_{pl}$ applied on the lungs surface. This pressure varies under the action of the diaphragm from -5 cmH₂O (i.e. -0.5 kPa) at end-exhalation to -8 cmH₂O (i.e. -0.8 kPa) at standard end-inhalation, and can reach -30 cmH₂O (i.e. -3 kPa) at full inhalation. The associated boundary traction on $\gamma$ is $\underline{t} = -p_{pl}\underline{n}$, where $\underline{n}$ denotes the outward unit normal vector. It corresponds to the pulled back boundary traction $\underline{T} = -p_{pl}J\underline{F}^{-T}\cdot\underline{N}_0$ on $\Gamma_0$, with $\underline{N}_0$ denoting the unit outward normal vector on $\Gamma_0$.

(iii) Interactions between the rib cage and the lungs are modeled by a bilateral contact between the lungs external surface and the rib cage surface. This contact is assumed to be frictionless, taking into account the lubrication by the pleural liquid. As the contact is bilateral, it is permanent, which means that surface points that are in contact initially remain so throughout breathing. These interactions are represented by a boundary traction, denoted by $\underline{T}_{c,rb}$, applied on the lung surface in contact with the rib cage $\Gamma_{0,c,rb}$. Since contact is frictionless, this traction is also associated with a pressure, i.e. $\underline{T}_{c,rb} = -p_{c,rb}J\underline{F}^{-T}\cdot\underline{N}_0$.

(iv) Rib cage displacements induced by intercostal muscles during breathing are different for everyone. They are imposed in the model through a prescribed motion of the rib cage surface. Patient-specific motion can be extracted from biomedical images, or generic motion (including no motion, i.e. fixed rib cage surface) can be used.

(v) Gravity loading, denoted $g$, can be applied considering any position, e.g. supine, prone or standing/seated position.

The virtual work of external forces is then

$$ W_{\text{ext}}(\underline{U}, \underline{U}^*) = - \int_{\Gamma_{0,c,in}} p_{\text{pl}} J (\underline{F}^{-T} \cdot \underline{N}_0) \cdot \underline{U}^* \, d\Gamma_0 \\ - \int_{\Gamma_{0,c,rb}} p_{\text{c,rb}} J (\underline{F}^{-T} \cdot \underline{N}_0) \cdot \underline{U}^* \, d\Gamma_0 + \int_{\Omega_0} \bar{\rho} \underline{g} \cdot \underline{U}^* \, d\Omega_0. \qquad (38) $$

Note that since we are imposing the thorax displacements the pressure quantity $p_{c,rb}$ will not be prescribed, but has the nature of a reaction force associated with frictionless contact. Nevertheless it is meaningful to include it in the external virtual work expression, as in practice this type of contact will often be modeled by means of penalization or mixed formulation, in which the Lagrange multiplier would be this pressure itself.

**Mechanical ventilation** In the case of mechanical ventilation, the driving force is now internal, i.e. given by the ventilator pressure, and we need to model boundary conditions representing the fact that the surrounding structures resist lung inflation.

(i) Air pressure is controlled by the ventilator. We again assume that air pressure is homogeneous, i.e. $p_f = p_{\text{vent}}$ with $p_{\text{vent}}$ the ventilated pressure. The associated stresses appear in Eq. (22) as a part of the mixture stresses. Using the homogeneous pressure assumption, the corresponding
volume integral can be transformed into a surface expression by Green's formula, as follows

$$ \int_{\Omega_0} p_f J \underline{C}^{-1} : (d\underline{U} \underline{E} \cdot \underline{U}^*) d\Omega_0 = \int_{\gamma} p_f \underline{n} \cdot \underline{U}^* d\gamma = \int_{\Gamma_0} p_f J (\underline{F}^{-T} \cdot \underline{N}_0) \cdot \underline{U}^* d\Gamma_0. \quad (39) $$

It means that the fluid pressure, acting on the whole lung volume, can be equivalently considered – under the homogeneous pressure assumption – as a following pressure applied on the lung surface $\gamma$.

(ii) Similarly to the interaction with the rib cage, interaction between the internal organs and the lungs involve a bilateral contact between the lungs external surface and the inner surface. This contact is frictionless, as the pleural liquid acts as a lubricant, hence has an associated reaction pressure applied on the lung external surface in contact with the inner surface $\Gamma_{0,c,in}$, which we denote by $p_{c,in}$.

The virtual work of external forces is then

$$ W_{\text{ext}}(\underline{U}, \underline{U}^*) = - \int_{\Gamma_{0,c,in}} p_{c,in} J (\underline{F}^{-T} \cdot \underline{N}_0) \cdot \underline{U}^* d\Gamma_0 \\ - \int_{\Gamma_{0,c,rb}} p_{c,rb} J (\underline{F}^{-T} \cdot \underline{N}_0) \cdot \underline{U}^* d\Gamma_0 + \int_{\Omega_0} \bar{\rho} \underline{g} \cdot \underline{U}^* d\Omega_0, \quad (40) $$

where $p_{c,in}$ should be modeled by Robin-type boundary conditions to represent surrounding organs that oppose lung inflation, with a stiffness coefficient to be adjusted according to the stiffness of each structure considered. Note that in this case $p_{c,rb}$ can be modeled in a similar manner with a large stiffness coefficient associated with the rib cage, since the rib cage motion is now in reaction to lung inflation, unlike for spontaneous breathing where it is a driving factor thereof. The driving factor is now, again, the internal fluid pressure, which is not included in $W_{\text{ext}}$ as it is already part of the mixture stress tensor $\underline{\Sigma}$, hence of the virtual work of internal forces $W_{\text{int}}$.

**Remarks on pleural pressure** Several studies have shown that pleural pressure is not homogeneous in the whole pleural space. For instance, in [Kallet, 2015] a pleural pressure gradient of about 7.5 cmH₂O has been reported between the apex and the base in upright position i.e. the pleural pressure is 7.5 cmH₂O more negative at the apex than the base.

The main explanation generally given for this gradient is gravity [Agostoni, 1972; Millar and Denison, 1989], which induces a smaller pressure, i.e. more negative, at the apex (which is at the top) than at the base (which is at the bottom) in standing position. If the pleural liquid was in an hydrostatic equilibrium, the gradient of pleural pressure would be 1 cmH₂O/cm [Agostoni, 1972; Lai-Fook, 2004], i.e. a gradient of 20 cmH₂O for a 20 cm high lung between the apex and the base.

The difference between the observed gradient and the hydrostatic gradient implies the existence of an intrinsic gradient that reduces the gravity effect. It can possibly be explained by the existence of contact areas between pleural membranes due to the presence of fatty tissue layers on the ribs, between ribs and the parietal pleura. These fatty tissue layers, 1 or 2 mm thick [Im et al., 1989], induce a tighter contact between parietal and visceral pleura, which may disrupt the hydrostatic equilibrium where barriers are formed in the pleural space [Lai-Fook, 2004]. This explanation is consistent with the fact that the gradient of pleural pressure between apex and base is also present in supine position. Since the main driver of lung inflation lies in the traction of the diaphragm,
it directly dictates the negative pleural pressure in the basal area, but this pressure is only partly propagated to the other parts of the lung surface.

The role of the pleural liquid circulation in the pleural pressure gradient has also been been raised. Indeed, the pleural liquid is highly viscous and is always circulating in the pleural space, in particular as part of a renewal mechanism. The pleural liquid is produced by capillary filtering in the upper parietal pleura and is resorbed by lymphangions in the basal, posterior costal parietal pleura and the diaphragmatic and mediastinal parietal pleura, inducing a downward pleural liquid flow that could be associated with a pressure drop along the flow [Wang, 1975; Miserocchi, 1997; Zocchi, 2002]. However, the complete turnover of pleural liquid is performed only twice a day, since about 15 mL are produced each day and the volume of pleural liquid is about 8 mL [Noppen et al., 2000], and therefore the pleural cavity can be reasonably considered as closed at the time scale of the breathing cycle, and therefore this explanation appears to be very unlikely.

In our modeling, the pleural pressure will be first assumed to be homogeneous on the whole lung surface. Gradients of pleural pressure will also be considered to take into account the non-homogeneity of the pleural pressure: from apex to base in standing position, from back to stomach in supine or prone position.

## 2.4 End-exhalation stress field and unloaded configuration

Since in general the initial geometry of the lungs (which may be obtained from biomedical images or from a generic model) is not unloaded as was already mentioned in Section 2.2, the unloaded configuration needs to be computed in order to take into account the initial stress field present in the initial configuration. This is a classical problem in biomechanics, since many living tissue contain prestresses due to a loaded physiological configuration [Gee et al., 2010; Sellier, 2011; Rausch et al., 2017] or inherent residual stresses [Fung, 1981; Genet et al., 2015b; Genet, 2019].

Let us consider that a pleural pressure $p_{pl,I}$ is applied on the configuration $\omega_I$ constructed from an image $I$. We aim to find the unloaded configuration $\Omega_0$ so that $\Omega_0$ is transformed into $\omega_I$ when $p_{pl,I}$ is applied on it. For the sake of simplicity, this unloaded configuration will serve as reference configuration for all subsequent computations. The boundary conditions to describe the transformation from the unloaded configuration $\Omega_0$ to a real configuration $\omega_I$ are simplified compared to the boundary conditions used to describe breathing, since only a negative pressure applied on the whole lung surface, which stands for the pleural pressure, is considered (see Figure 4). In addition to the negative pleural pressure, the rigid body degrees of freedom are blocked through the following procedure: the four nodes (one origin node & three directional nodes) defining the most orthogonal trihedron are found within the mesh, the origin is blocked in all three directions, the first directional node is blocked in the second & third directions, and the second directional node is blocked in the third direction.

In this problem, a configuration $\omega_I$ is given and the unknowns are the transformation $\chi^{-1}(\underline{x})$ (or equivalently the displacement $\underline{u}(\underline{x})$) from $\omega_I$ to the unloaded configuration $\Omega_0$, as well as the reference porosity on the deformed configuration $\phi_{f0}(\underline{x}) = \frac{dV_{f0}}{dV}$. The boundary conditions applied from $\Omega_0$ to $\omega_I$ are given. The problem is formulated following [Govindjee and Mihalic, 1998], where a change of variable is performed from $\underline{X}$ to $\underline{x}$, which is extended to poromechanics for the first time in this paper. We denote with small letters the quantities defined on the known deformed
<figure><img src="image_7.png" /><figcaption>Figure 4: Boundary conditions used to compute the unloaded configuration $\Omega_0$, then the deformed configuration $\omega_i$ from the initial configuration $\omega_e$. In addition to preventing rigid body motion, only pleural pressure $p_{pl,e}$ is considered for the transformation $\chi_0$ between the unloaded and the end-exhalation configurations. However, the breathing transformation $\chi_b$ between the end-exhalation to the end-inhalation configurations uses more complex boundary conditions, which depend on the type of breathing. In the case of spontaneous breathing, we consider the rib cage displacements $\underline{U}_{thorax}$, the pleural pressure at end-inhalation $p_{pl,i}$, a frictionless bilateral contact between the rib cage surface and lungs surface and the gravity volume loading. In the case of mechanical ventilation, an inner surface subjected to Robin's conditions is also in contact with the inner lung surface.</figcaption></figure>
configuration $\omega_I$ (see Figure 4), such as:

$$ \left\{ \begin{array}{l} \underline{u} := -\underline{U} \circ \underline{\chi}^{-1} \\ \underline{f} := \underline{F}^{-1} \circ \underline{\chi}^{-1} \\ j := J^{-1} \circ \underline{\chi}^{-1} \\ \phi_{f0} := (J^{-1} \circ \underline{\chi}^{-1}) (\Phi_{f0} \circ \underline{\chi}^{-1}) \end{array} \right. \qquad (41) $$

We now use the principle of virtual work written in the deformed configuration, i.e.

$$ \forall \underline{u}^*, \int_{\omega_I} \underline{\sigma}(\underline{u}) := \underline{\varepsilon}(\underline{u}^*) d\omega = - \int_{\gamma_I} p_{\text{pl},I} \underline{n} \cdot \underline{u}^* d\gamma, \qquad (42) $$

where the Cauchy stress tensor $\underline{\sigma}$ is written after changing variables as

$$ \begin{aligned} \underline{\sigma}(\underline{u}(\underline{x})) &= \frac{1}{J(\underline{\chi}^{-1}(\underline{x}))} \underline{F}(\underline{\chi}^{-1}(\underline{x})) \cdot \underline{\Sigma}(\underline{\chi}^{-1}(\underline{x})) \cdot \underline{F}^T(\underline{\chi}^{-1}(\underline{x})) \\ &= j \underline{f}^{-1}(\underline{x}) \cdot \underline{\Sigma}(\underline{\chi}^{-1}(\underline{x})) \cdot \underline{f}^{-T}(\underline{x}) \end{aligned} \qquad (43) $$

and the linearized strain tensor $\underline{\varepsilon}$ is

$$ \underline{\varepsilon}(\underline{u}) := (\underline{\nabla} \underline{u})_{\text{sym}} = \frac{1}{2} (\underline{\nabla} \underline{u} + \underline{\nabla} \underline{u}^T). \qquad (44) $$

Note that it is convenient here to write the principle of work in the deformed configuration since it is known, instead of the reference configuration as in (14). This is why it is the Cauchy stress tensor and the linearized strain tensor that appear, instead of the second Piola-Kirchhoff stress tensor and Green-Lagrange strain tensor. The expressions of the second Piola-Kirchhoff stress tensor and the fluid pressure can then be derived by introducing the same change of variables into the general expressions (22) & (23).

In this paper, we choose to estimate the unloaded configuration from the end-exhalation configuration, as shown in Figure 4.

## 2.5 Summary of main equations

### 2.5.1 General poromechanics formulation for lung modeling

The problem to be solved to determine the deformation of lungs in our poromechanical framework is described by the following system

$$ \text{Find } (\underline{U}, \underline{\Phi}_f) / \left\{ \begin{array}{l} \forall \underline{U}^*, \int_{\Omega_0} \frac{\partial \bar{\psi}}{\partial \underline{E}} : (d\underline{U} \underline{E} \cdot \underline{U}^*) d\Omega_0 = W_{\text{ext}}(\underline{U}, \underline{U}^*) \\ p_f = - \frac{\partial \bar{\psi}_s}{\partial \underline{\Phi}_s} \end{array} \right. \qquad (45) $$

In this system, the two equations are coupled with the two unknowns $\underline{U}$ (i.e. the mixture displacement) and $\underline{\Phi}_f$ (i.e. the deformed porosity by unit reference mixture volume). The first equation
governs the global mechanical equilibrium, and the second one expresses the local hydrostatic pressure equilibrium.

Using the decomposition of $\frac{\partial \bar{\psi}}{\partial \underline{E}}$ given in Eq. (22) as well as the transformation of the fluid pressure term with the homogeneous fluid pressure assumption, the first equation of the system gives

$$ \text{Find } (\underline{U}, \Phi_f) / \left\{ \begin{array}{l} \forall \underline{U}^*, \int_{\Omega_0} \frac{\partial \bar{\psi}_s}{\partial \underline{E}} : (d_{\underline{U}} \underline{E} \cdot \underline{U}^*) d\Omega_0 = W_{\text{ext}}(\underline{U}, \underline{U}^*) + \int_{\Gamma_0} p_f J(\underline{F}^{-T} \cdot \underline{n}) \cdot \underline{U}^* d\Gamma_0 \\ p_f = -\frac{\partial \bar{\psi}_s}{\partial \Phi_s} \end{array} \right. \quad (46) $$

It is interesting to notice that the internal fluid pressure has the same effect as a following pressure applied on the boundary of the mixture.

With the choice of $\bar{\psi}_s$ that was made in Eq. (24), we have

$$ \text{Find } (\underline{U}, \Phi_f) / \left\{ \begin{array}{l} \forall \underline{U}^*, \int_{\Omega_0} \frac{\partial \bar{W}_{\text{skel}}(\underline{E})}{\partial \underline{E}} : (d_{\underline{U}} \underline{E} \cdot \underline{U}^*) d\Omega_0 = W_{\text{ext}}(\underline{U}, \underline{U}^*) + \int_{\Gamma_0} p_f J(\underline{F}^{-T} \cdot \underline{n}) \cdot \underline{U}^* d\Gamma_0 \\ p_f = -\frac{\partial \bar{W}_{\text{bulk}}}{\partial \Phi_s} \end{array} \right. \quad (47) $$

and we see that the coupling is weak, since the first equation – similar to that of a hyperelastic solid – is independent of the porosity $\Phi_f$, which can be computed in a post-processing step using the second equation.

With the particular choice made in Eq. (29) for $\bar{W}_{\text{bulk}}$, we even have a closed form expression for the porosity

$$ \Phi_f = J - \Phi_s \quad \text{with} \quad \Phi_s = \frac{1}{\frac{p_f}{\bar{\kappa}} + \frac{1}{\Phi_{s0}}}, \quad (48) $$

with $\Phi_{s0} := 1 - \Phi_{f0}$. In the case of free breathing where $p_f = 0$ – or when we tend to an incompressible behavior for the solid material, i.e. $\bar{\kappa}$ is very large – the quantity $\Phi_s$ is constant and equal to $\Phi_{s0}$, which leads to the geometric identity

$$ 1 - \Phi_{f0} = J \left(1 - \phi_f \circ \underline{\chi}\right). \quad (49) $$

### 2.5.2 Formulation with the additional energy $\bar{W}_{\text{por}}$

When using an additional energy $\bar{W}_{\text{por}}$, recalling (32)–(33) the formulation is then

$$ \text{Find } (\underline{U}, \Phi_f) / \left\{ \begin{array}{l} \forall \underline{U}^*, \int_{\Omega_0} \left( \frac{\partial \bar{W}_{\text{skel}}(\underline{E})}{\partial \underline{E}} + \frac{\partial \bar{W}_{\text{bulk}}(\Phi_s)}{\partial \Phi_s} J \underline{C}^{-1} \right) : (d_{\underline{U}} \underline{E} \cdot \underline{U}^*) d\Omega_0 = W_{\text{ext}}(\underline{U}, \underline{U}^*) \\ p_f = -\frac{\partial \bar{W}_{\text{bulk}}(\Phi_s)}{\partial \Phi_s} + \frac{\partial \bar{W}_{\text{por}}(\Phi_f)}{\partial \Phi_f} \end{array} \right. \quad (50) $$
The coupling is strong here, since the first equation depends on the current porosity $\Phi_f = J - \Phi_s$. However, this coupling is only local since the second equation – where $p_f$ is known – provides a relationship between $J$ and $\Phi_f$.

### 2.5.3 Contact-like formulation

Recalling Eqs. (34)–(37), the contact-like formulation to ensure a positive porosity – formulated using the contact pressure $p_c$ as an extra unknown – can be summarized as follows

$$ \text{Find } (\underline{U}, \Phi_f, p_c) / \left\{ \begin{array}{l} \forall \underline{U}^*, \int_{\Omega_0} \left( \frac{\partial \bar{\psi}_s}{\partial \underline{E}} - p J \underline{C}^{-1} \right) : (d\underline{U} \underline{E} \cdot \underline{U}^*) d\Omega_0 = W_{\text{ext}}(\underline{U}, \underline{U}^*) \\ p = p_f \mathcal{H}(\Phi_f) + p_c (1 - \mathcal{H}(\Phi_f)) = - \frac{\partial \bar{\psi}_s}{\partial \Phi_s} \\ \Phi_f \cdot p_c = 0 \\ \Phi_f \ge 0 \\ p_c \ge 0 \end{array} \right. \qquad (51) $$

A more compact form can be found by introducing the auxiliary unknown $\hat{\Phi}_f$ such that

$$ p_f = - \frac{\partial \bar{\psi}_s}{\partial \Phi_s} (\underline{E}, \Phi_s = J - \hat{\Phi}_f). $$

This quantity coincides with $\Phi_f$ when $\hat{\Phi}_f \ge 0$, hence

$$ \Phi_f = \hat{\Phi}_f \mathcal{H}(\hat{\Phi}_f), \qquad (52) $$

and therefore we have the following equivalent formulation

$$ \text{Find } (\underline{U}, \hat{\Phi}_f) / \left\{ \begin{array}{l} \forall \underline{U}^*, \int_{\Omega_0} \left( \frac{\partial \bar{\psi}_s}{\partial \underline{E}} + \frac{\partial \bar{\psi}_s}{\partial \Phi_s}(\underline{E}, J - \hat{\Phi}_f \mathcal{H}(\hat{\Phi}_f)) J \underline{C}^{-1} \right) : (d\underline{U} \underline{E} \cdot \underline{U}^*) d\Omega_0 = W_{\text{ext}}(\underline{U}, \underline{U}^*) \\ p_f = - \frac{\partial \bar{\psi}_s}{\partial \Phi_s}(\underline{E}, J - \hat{\Phi}_f) \end{array} \right. \qquad (53) $$

This problem has a similar structure to that of (50), i.e. $J$ and $\hat{\Phi}_f$ are locally coupled by the second equation. However, here $\hat{\Phi}_f$ is only used in the first equation when it is positive, i.e. when $\Phi_f = \hat{\Phi}_f$.

As is well-known with contact formulations, in practice the convergence of associated numerical procedures can be difficult to achieve. Typically here, difficulties might arise due to the very different stiffness levels between parts where the internal contact is activated or not, especially when the solid material is rather incompressible.

### 2.5.4 Inverse problem

In the case of the inverse problem to estimate the unloaded configuration, a loaded configuration $\omega_I$ – which can be taken from an image - is given, as well as the associated deformed porosity $\phi_f$.
The two unknowns are the inverse displacement on the deformed configuration $\underline{u}$ and the reference porosity on the deformed configuration $\phi_{f0}$.

In the case of the penalized strategy presented in Section 2.1.4, we have the following formulation:

Find $(\underline{u}, \phi_{f0}) /$

$$
\left\{
\begin{array}{l}
\forall \underline{u}^*, \quad \displaystyle\int_{\omega_1} \underline{\sigma}(\underline{u}, \phi_{f0}) : \underline{\varepsilon}(\underline{u}^*) d\omega = - \displaystyle\int_{\gamma_1} p_{pl} \underline{n} \cdot \underline{u}^* d\gamma \\
\underline{\sigma}(\underline{u}, \phi_{f0}) = j \underline{f}^{-1} \cdot \underline{\Sigma}(\underline{u}, \phi_{f0}) \circ \underline{\chi}^{-1} \cdot \underline{f}^{-T} \\
\underline{\Sigma}(\underline{u}, \phi_{f0}) = \frac{\partial \bar{\psi}_s}{\partial \underline{E}} (\underline{E}, \Phi_s = (1-\phi_f) \circ \underline{\chi} \cdot J; \Phi_{f0} = \phi_{f0} \circ \underline{\chi} \cdot J) \qquad (54) \\
\qquad + \frac{\partial \bar{\psi}_s}{\partial \Phi_s} (\underline{E}, \Phi_s = (1-\phi_f) \circ \underline{\chi} \cdot J; \Phi_{f0} = \phi_{f0} \circ \underline{\chi} \cdot J) J \underline{C}^{-1} \\
p_f = - \frac{\partial \bar{\psi}_s}{\partial \Phi_s} (\underline{E}, \Phi_s = (1-\phi_f) \circ \underline{\chi} \cdot J; \Phi_{f0} = \phi_{f0} \circ \underline{\chi} \cdot J)
\end{array}
\right.
$$

Note that the $\Phi_{f0}$ dependency of $\bar{\psi}_s$ was not explicitly written before, since in the direct problem $\Phi_{f0}$ is a given parameter, not a variable as in the inverse problem. Nevertheless, positive porosity in the unloaded configuration will be automatically enforced due to the specifications prescribed on the penalization term.

In order to ensure a positive reference porosity using the contact-like formulation introduced in Section 2.1.4, we introduce the auxiliary unknown $\hat{\phi}_{f0}$ such that

$$ p_f = - \frac{\partial \bar{\psi}_s}{\partial \Phi_s} (\underline{E}, \Phi_s = (1-\phi_f) \circ \underline{\chi} \cdot J; \Phi_{f0} = \hat{\phi}_{f0} \circ \underline{\chi} \cdot J). \quad (55) $$

Note that as for the previous formulation, the $\Phi_{f0}$ dependency is made explicit as for inverse problems it is a variable, not a parameter, of the problem. This auxiliary variable can be negative; however, it coincides with the actual reference porosity $\phi_{f0}$ when it is positive, i.e.

$$ \phi_{f0} = \hat{\phi}_{f0} \mathcal{H}(\hat{\phi}_{f0}). \qquad (56) $$

However, the fluid pressure should satisfy (23) with the actual value of $\Phi_{f0} = \phi_{f0} \circ \underline{\chi} \cdot J$, which is in contradiction with (55). Therefore we need to correct the solid free energy, which we propose to do by considering the modified solid free energy

$$ \bar{\psi}_s^c = \bar{\psi}_s - C_p(\Phi_s - J), \qquad (57) $$

with

$$ C_p = - \frac{\partial \bar{\psi}_s}{\partial \Phi_s} (\underline{E}, (1-\phi_f) \circ \underline{\chi} \cdot J; \hat{\phi}_{f0} \circ \underline{\chi} \cdot J) + \frac{\partial \bar{\psi}_s}{\partial \Phi_s} (\underline{E}, (1-\phi_f) \circ \underline{\chi} \cdot J; \phi_{f0} \circ \underline{\chi} \cdot J), $$

a pressure correction designed so that

$$
\begin{align*}
p_f &= - \frac{\partial \bar{\psi}_s^c}{\partial \Phi_s} (\underline{E}, (1-\phi_f) \circ \underline{\chi} \cdot J; \phi_{f0} \circ \underline{\chi} \cdot J) \\
&= - \frac{\partial \bar{\psi}_s}{\partial \Phi_s} (\underline{E}, (1-\phi_f) \circ \underline{\chi} \cdot J; \phi_{f0} \circ \underline{\chi} \cdot J) + C_p \\
&= - \frac{\partial \bar{\psi}_s}{\partial \Phi_s} (\underline{E}, (1-\phi_f) \circ \underline{\chi} \cdot J; \hat{\phi}_{f0} \circ \underline{\chi} \cdot J).
\end{align*}
$$
Note that we expect $-\frac{\partial \bar{\psi}_s}{\partial \Phi_s}(\underline{E}, \Phi_s; \Phi_{f0})$ to be increasing when $\Phi_{f0}$ decreases, by construction of the pressure law (23). In other words, it should require more fluid pressure to compress the solid material to a given fraction $\Phi_s$ in the deformed configuration when the solid fraction in the reference configuration is larger – or equivalently when the fluid fraction $\Phi_{f0}$ is smaller. Therefore the pressure correction $C_p$ should be positive when it is activated, since this occurs when $\Phi_{f0} = 0$ and $\hat{\phi}_{f0} < 0$. Furthermore, the correction applied in (57) is designed to have no incidence on the stress tensor, and we therefore have the following formulation

$$ \text{Find } (\underline{u}, \hat{\phi}_{f0}) / \left\{ \begin{array}{l} \forall \underline{u}^*, \int_{\omega_1} \underline{\sigma}(\underline{u}, \hat{\phi}_{f0}) : \underline{\varepsilon}(\underline{u}^*) d\omega = - \int_{\gamma_1} p_{pl} \underline{n} \cdot \underline{u}^* d\gamma \\ \underline{\sigma}(\underline{u}, \hat{\phi}_{f0}) = j \underline{f}^{-1} \cdot \underline{\Sigma}(\underline{u}, \hat{\phi}_{f0}) \circ \underline{\chi}^{-1} \cdot \underline{f}^{-T} \\ \underline{\Sigma}(\underline{u}, \hat{\phi}_{f0}) = \frac{\partial \bar{\psi}_s}{\partial \underline{E}}(\underline{E}, \Phi_s = (1-\phi_f) \circ \underline{\chi} \cdot J; \Phi_{f0} = \hat{\phi}_{f0} \mathcal{H}(\hat{\phi}_{f0}) \circ \underline{\chi} \cdot J) \\ \qquad + \frac{\partial \bar{\psi}_s}{\partial \Phi_s}(\underline{E}, \Phi_s = (1-\phi_f) \circ \underline{\chi} \cdot J; \Phi_{f0} = \hat{\phi}_{f0} \mathcal{H}(\hat{\phi}_{f0}) \circ \underline{\chi} \cdot J) J \underline{C}^{-1} \\ p_f = - \frac{\partial \bar{\psi}_s}{\partial \Phi_s}(\underline{E}, \Phi_s = (1-\phi_f) \circ \underline{\chi} \cdot J; \Phi_{f0} = \hat{\phi}_{f0} \circ \underline{\chi} \cdot J) \end{array} \right. \quad (58) $$

This problem is strongly coupled since $\hat{\phi}_{f0}$ – which is determined by the last equation – is also present in the first equation. However, as for the direct problem, this coupling is only local. In addition, note that the above correction in the solid free energy is rather ad hoc – since it depends on the solution of the inverse problem itself – albeit in practice, in the pulmonary setting, the extent of the regions in which the correction is activated is expected to be extremely restricted, as even in the unloaded state the porosity is not expected to be close to zero in significantly large regions.

## 2.6 Implementation aspects

### 2.6.1 Geometric model and unloaded configuration

The rib cage and inner surfaces generation uses the lung geometries. First, the lung surface is used to compute a binary image. Then, for each surface, a second binary image is performed by drawing contours manually with the software MeVisLab in order to get an extension of the lungs. Finally, surfaces are extracted from each of the binary images and meshed with the software Gmsh to get the rib cage surface and the inner surface.

In order to compute the unloaded configuration, the material model and method described in Section 2.4 have been implemented in Python using the FEniCS library [Alnæs et al., 2015; Logg et al., 2012], and the resulting code is freely available<fnref n="1" />. A custom Newton-Raphson procedure is used to solve the nonlinear problem.

<footnote n="1">¹https://gitlab.inria.fr/mgenet/dolfin_mech</footnote>
### 2.6.2 Finite element solution for deformed configuration

The deformed configuration is computed with the Abaqus software. We now present the details of the model implementation.

**Specific finite elements** The lungs are meshed with linear tetrahedral volume elements, whereas both the rib cage surface and the inner surface are meshed with linear triangle surface elements. The corresponding finite elements are respectively called C3D4 and S3 in Abaqus.

**Mechanical properties** The skeleton energy $\overline{W}_{skel}$ is written in an Abaqus subroutine UHYPER. The local equilibrium of the mixture depending on the bulk energy $\overline{W}_{bulk}$ is also implemented in the subroutine in order to compute the deformed porosity as a state variable. The inner surface is defined with a Neo-Hookean material, with shear modulus $2 \cdot 10^3$ kPa and bulk modulus 6.67 kPa.

**Contact boundary conditions** Concerning the contact with the thorax, a master-slave, finite-sliding, node-to-surface contact pair formulation is used, with the rib cage as master and the lung surface as slave. The surfaces are adjusted at the start of the simulation in order to remove gaps and overlaps. The contact is defined as a frictionless, hard bilateral contact which does not allow any separation of surfaces once it is established. In the case of the mechanical ventilation, the contact with the inner surface has the same properties as the contact with the thorax except that separation between surfaces is allowed.

**Robin boundary conditions** Robin boundary conditions are implemented on the inner surface using spring elements. Each node of the inner surface is linked to three linear spring elements – one for each displacement component – of type SPRING1, attaching the node to the end-exhalation position of that node. The spring stiffness $k_{spring}$ is different for each node in order to take into account the size of the surface elements connected to this node. The spring stiffness is then $k_{spring} = \lambda_{node}K$, where $\lambda_{node}$ is the coefficient of the lumped mass matrix associated with the node, and $K$ is a reference stiffness coefficient.

**Solution process** The deformed configuration is computed in two steps. First, the stresses and the strains in the initial configuration are computed using the reference nodes coordinates and the reference porosity of the unloaded configuration as computed before. It uses the same model as for the estimation of the unloaded configuration. It is checked that the nodes coordinates are the same as in the initial configuration from which the unloaded configuration was estimated. Then, the deformed configuration is computed with the more complex boundary conditions, including contact. Both problems are solved using the Newton-Raphson procedure built into Abaqus.

## 3 Illustrative results

We now present some illustrative results using our pulmonary poromechanical model, with test cases motivated by clinical applications and patient-specific modeling. The end-exhalation pleural pressure $p_{pl,e}$ and the end-inhalation pleural pressure $p_{pl,i}$ are chosen in the physiological range, and the mechanical parameters are chosen in order to reproduce the volumetric response to a change of pressure as observed in experimental data [Gibson and Pride, 1976] (see Figure 2). In all cases,
the thorax is fixed, and gravity is neglected. For each illustrative case, the parameters used are reported in the corresponding table.

## 3.1 Geometries

All the computations in this section use a generic right lung mesh obtained based on the Zygote geometric model<fnref n="2" /> in the end-exhalation configuration, and a thorax mesh generated from the lung mesh in order to have the thorax surface and the lung surface in contact in the initial configuration. The thorax corresponds to the hard structures surrounding lungs, mainly composed of the ribs and the spine. The lung mesh contains 4052 tetrahedral elements, whereas the thorax and the inner surface comprise 2678 and 2260 triangular elements respectively. Both lung and thorax meshes are shown in Figure 5.

<figure><img src="image_8.png" /><figcaption>Figure 5: Lung and thorax meshes in the initial configuration, i.e. end-exhalation. (Left) Lung mesh obtained based on Zygote model<fnref n="2" />. (Middle) Two views of the thorax mesh, lateral view (top) and top view (bottom). (Right) Both meshes together in a coronal view. The lung surface is in contact with the thorax surface.</figcaption></figure>

<footnote n="2">²http://www.zygote.com</footnote>
## 3.2 Influence of the positive porosity constraint on the unloaded configuration

<table><caption>Table 1: Parameters used for the breathing simulations when studying the influence of the porosity constraint on the unloaded configuration. Top rows correspond to the loading parameters, i.e. fluid pressure and pleural pressures, whereas bottom rows give the coefficients of the constitutive behavior of the pulmonary mixture. The solid bulk modulus is not needed here since $p_f = 0$.</caption><thead><tr><th>Name</th><th>Symbol</th><th>Unit</th><th>Value</th></tr></thead><tbody><tr><td>Air pressure</td><td>p<sub>f</sub></td><td>[kPa]</td><td>0</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 rowspan="4">Skeleton energy parameters</td><td>α̃</td><td>[kPa]</td><td>1.6 · 10<sup>-1</sup></td></tr><tr><td>δ</td><td>[-]</td><td>0.5</td></tr><tr><td>β̃<sub>1</sub></td><td>[kPa]</td><td>0.2</td></tr><tr><td>β̃<sub>2</sub></td><td>[kPa]</td><td>0.4</td></tr><tr><td>Solid bulk modulus</td><td>κ̄</td><td>[kPa]</td><td>-</td></tr></tbody></table>

The first step to simulate breathing is to estimate the unloaded configuration. As explained in Sections 2.1.4 and 2.4, the estimation of the unloaded configuration is associated with a problem concerning the porosity to be kept positive. We study in this paragraph the effect of using the contact-like formulation in a realistic case, which is the asymptotic limit of the pore energy formulation.

An initial end-exhalation porosity distribution $\phi_{f,e}$ was generated according to a Gaussian distribution (mean: 0.5, std: 0.13) as shown in grey in Figure 6. The porosity values are randomly distributed in space over the lung volume. The pressure loading and material parameters used in the simulations are given in Table 1.

Two simulations are presented: one without any constraint on the porosity using the model of Equation (24), the other using the contact-like formulation proposed in Equation (34). The resulting reference porosity $\Phi_f$ in both cases is then computed as displayed in Figure 6, in blue and red, respectively.

## 3.3 Influence of the stiffness

We present here an illustration of the model in the free breathing case, and the fluid pressure is then zero. The unloaded configuration is first computed from the initial configuration using the inverse problem described in Section 2.4. Then, the deformed configuration is computed using the complex boundary conditions described in Section 2.3: a pleural pressure is applied on the lung surface and makes the lung inflate and slide against the thorax surface; the thorax displacement is prescribed (here it is null). The loading and material parameters used for the computations are provided in Table 2. The same material parameters are used for both the estimation of the unloaded configuration and the computation of the deformed configuration.

The impact of the stiffness parameter $\bar{\alpha}$, which is part of the skeleton energy presented in Equation (27), is studied using three different values: $8 \cdot 10^{-1}$ kPa, $8 \cdot 10^{-2}$ kPa, $8 \cdot 10^{-3}$ kPa. The
<figure><img src="image_9.png" /><figcaption>Figure 6: Comparison of the reference porosity field computed without or with the contact-like formulation. (Left) Given end-exhalation porosity shown on the end-exhalation mesh. (Right) Distribution of the given end-exhalation porosity plotted in grey in the histogram. The red and blue distributions are the computed reference porosity with and without the contact-like formulation, respectively. When not using any strategy to ensure a positive porosity, the computed reference porosity of some cells is negative, whereas the computed reference porosity remains positive or zero when using the contact-like formulation.</figcaption></figure>

<table><caption>Table 2: Parameters used for the breathing simulations when studying the stiffness impact. Top rows correspond to the loading parameters, fluid pressure and pleural pressures, whereas bottom rows feature the coefficients of the constitutive behavior of the pulmonary mixture. The solid bulk modulus is not needed since $p_f = 0$.</caption><thead><tr><th>Name</th><th>Symbol</th><th>Unit</th><th>Value</th></tr></thead><tbody><tr><td>Air pressure</td><td>p<sub>f</sub></td><td>[kPa]</td><td>0</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</td></tr><tr><td><strong>Skeleton energy parameters</strong></td><td></td><td></td><td></td></tr><tr><td>&alpha;</td><td>[kPa]</td><td>[8 &middot; 10<sup>-3</sup>, 8 &middot; 10<sup>-2</sup>, 8 &middot; 10<sup>-1</sup>]</td><td></td></tr><tr><td>&delta;</td><td>[-]</td><td></td><td>0.5</td></tr><tr><td>&beta;<sub>1</sub></td><td>[kPa]</td><td></td><td>0.1</td></tr><tr><td>&beta;<sub>2</sub></td><td>[kPa]</td><td></td><td>0.2</td></tr><tr><td>Solid bulk modulus</td><td>&kappa;</td><td>[kPa]</td><td>-</td></tr></tbody></table>
results are shown in Figure 7.

<figure><img src="image_10.png" /><figcaption>Figure 7: Volume response of lung to a change of pleural pressure in the free breathing case for different values of the stiffness parameter $\bar{\alpha}$ in kPa. (a) Lung volume change $J$ normalized by its value at end-exhalation, i.e. at the pressure 0.5 kPa, with respect to the pleural pressure loading. (b) The unloaded configurations (top row), are computed from the initial configurations (middle row) with simplified boundary conditions using the inverse problem presented in Section 2.4. Then, the deformed configurations (bottom row) are computed considering contact between the lungs and the thorax as described in Section 2.3.</figcaption></figure>
## 3.4 Comparison between free breathing and ventilated breathing

Mechanical ventilation can produce an injury called Ventilator-Induced Lung Injury (VILI) at the alveolar scale as a consequence of excessive pressure, which can induce over-expansion and large stresses in the alveoli. Clinicians are interested in improving mechanical ventilation to reduce such risks. To that aim they can tune several parameters like the Positive End-Expiratory Pressure (PEEP), the plateau pressure or the breathing rate. Lung models could be used to optimize these parameters, by allowing to quantify their impact on the parenchyma. For instance, our model allows to study the impact of fluid pressure induced by the ventilator on stresses, which we will now illustrate.

From Equation (25), we can decompose the mixture stress as follows:

$$ \underline{\underline{\Sigma}} = \underline{\underline{\Sigma}}_{\text{skel}} + \underline{\underline{\Sigma}}_{\text{p}} \quad \text{with} \quad \left\{ \begin{array}{l} \underline{\underline{\Sigma}}_{\text{skel}} := \frac{\partial \bar{W}_{\text{skel}}}{\partial \underline{\underline{E}}} \\ \underline{\underline{\Sigma}}_{\text{p}} := -p_{\text{f}} J \underline{\underline{C}}^{-1} \end{array} \right. \qquad (59) $$

where $\underline{\underline{\Sigma}}$ represents the total stress in the mixture (which equilibrates the boundary loadings on the mixture as shown by Equation (14)), $\underline{\underline{\Sigma}}_{\text{skel}}$ represents the skeleton stress (a function of the global deformation of the mixture that equilibrates the total pressure acting on the system, i.e. $p_{\text{tot}} = p_{\text{pl}} - p_{\text{f}}$, as shown by Equation (47)), and $\underline{\underline{\Sigma}}_{\text{p}}$ is the total hydrostatic stress. This mixture hydrostatic stress can be decomposed into solid and fluid parts, based on the porosity:

$$ \underline{\underline{\Sigma}}_{\text{p}} = \underline{\underline{\Sigma}}_{\text{p,s}} + \underline{\underline{\Sigma}}_{\text{p,f}} \quad \text{with} \quad \left\{ \begin{array}{l} \underline{\underline{\Sigma}}_{\text{p,s}} := (1 - \Phi_{\text{f}}) \underline{\underline{\Sigma}}_{\text{p}} \\ \underline{\underline{\Sigma}}_{\text{p,f}} := \Phi_{\text{f}} \underline{\underline{\Sigma}}_{\text{p}} \end{array} \right. \qquad (60) $$

Thus, the mixture stress can be decomposed into solid and fluid parts:

$$ \underline{\underline{\Sigma}} = \underline{\underline{\Sigma}}_{\text{s}} + \underline{\underline{\Sigma}}_{\text{f}} \quad \text{with} \quad \left\{ \begin{array}{l} \underline{\underline{\Sigma}}_{\text{s}} := \underline{\underline{\Sigma}}_{\text{skel}} + (1 - \Phi_{\text{f}}) \underline{\underline{\Sigma}}_{\text{p}} \\ \underline{\underline{\Sigma}}_{\text{f}} := \Phi_{\text{f}} \underline{\underline{\Sigma}}_{\text{p}} \end{array} \right. \qquad (61) $$

since $\underline{\underline{\Sigma}}_{\text{skel}}$ is a solid stress by definition.

<table><caption>Table 3: Typical range for the variations of pressure and fluid volume during a normal inhalation in both types of breathing [Goligher et al., 2016].</caption><thead><tr><th></th><th></th><th>Spontaneous breathing</th><th>Mechanical ventilation</th></tr></thead><tbody><tr><td>Δp<sub>pl</sub></td><td>[kPa]</td><td>+0.3</td><td>-</td></tr><tr><td>Δp<sub>f</sub></td><td>[kPa]</td><td>-</td><td>+3</td></tr><tr><td>ΔV<sub>air</sub></td><td>[L]</td><td>+0.5</td><td>+0.3</td></tr></tbody></table>

Both types of breathing – spontaneous breathing and mechanical ventilation – result in different stresses. Indeed, in the case of spontaneous breathing, the driving force of breathing is the variation of pleural pressure, whereas in the case of ventilated breathing it is the fluid pressure produced by the ventilator that acts to inflate both the lungs and the surrounding structures, including the rib
cage. The typical ranges of the pressure variation as well as the volume variation for both types of breathing are reported in Table 3. The variation of ventilation pressure is much larger than the variation of pleural pressure while the deformation is smaller, since the ventilation pressure must inflate the rib cage in addition to lungs, which also means that the pleural pressure will become positive, hence counterproductive. In the case of spontaneous breathing, we have

$$p_f^{\text{sb}} = 0, \qquad (62)$$

which entails

$$\left\{ \begin{array}{l} \Sigma_{\equiv f}^{\text{sb}} = \Sigma_{\equiv p}^{\text{sb}} = 0 \\ \Sigma_{\equiv s}^{\text{sb}} = \Sigma_{\equiv \text{skel}}^{\text{sb}} = \Sigma_{\equiv \text{skel}}^{\text{sb}}, \end{array} \right. \qquad (63)$$

whereas in the case of ventilated breathing we have

$$p_f^{\text{vb}} = p_{\text{vent}}, \qquad (64)$$

and therefore

$$\left\{ \begin{array}{l} \Sigma_{\equiv f}^{\text{vb}} = -\Phi_f p_{\text{vent}} J \underline{C}^{-1} \\ \Sigma_{\equiv s}^{\text{vb}} = \Sigma_{\equiv \text{skel}}^{\text{vb}} - (1 - \Phi_f) p_{\text{vent}} J \underline{C}^{-1}. \end{array} \right. \qquad (65)$$

The main difference lies in the contribution of the fluid pressure in the solid stress. Whereas it is zero in the case of spontaneous breathing, the fluid pressure is large in ventilated breathing. The skeleton stress $\Sigma_{\equiv \text{skel}}$, itself is slightly different as the deformation is somewhat larger in spontaneous breathing than in ventilated breathing.

<table><caption>Table 4: Parameters used for the comparison between spontaneous breathing and ventilated breathing. Top rows correspond to the loading parameters, fluid pressure and pleural pressures, whereas bottom rows feature the coefficients of the constitutive behavior of the pulmonary mixture. The solid bulk modulus is not needed since $p_f = 0$.</caption><thead><tr><th>Name</th><th>Symbol</th><th>Unit</th><th>Spontaneous breathing</th><th>Mechanical ventilation</th></tr></thead><tbody><tr><td><strong>End-inhalation air pressure</strong></td><td><em>p</em><sub>f</sub></td><td>[kPa]</td><td><strong>0</strong></td><td><strong>3</strong></td></tr><tr><td>End-exhalation pleural pressure</td><td><em>p</em><sub>pl,e</sub></td><td>[kPa]</td><td>-0.5</td><td>-0.5</td></tr><tr><td><strong>End-inhalation pleural pressure</strong></td><td><em>p</em><sub>pl,i</sub></td><td>[kPa]</td><td><strong>-0.8</strong></td><td>-</td></tr><tr><td rowspan="4">Skeleton energy parameters</td><td><em>α</em></td><td>[kPa]</td><td>8 &middot; 10<sup>-2</sup></td><td>8 &middot; 10<sup>-2</sup></td></tr><tr><td><em>δ</em></td><td>[-]</td><td>0.5</td><td>0.5</td></tr><tr><td><em>β</em><sub>1</sub></td><td>[kPa]</td><td>0.1</td><td>0.1</td></tr><tr><td><em>β</em><sub>1</sub></td><td>[kPa]</td><td>0.2</td><td>0.2</td></tr><tr><td>Solid bulk modulus</td><td><em>κ</em></td><td>[kPa]</td><td>-</td><td><strong>10<sup>2</sup></strong></td></tr></tbody></table>

The solid stress difference between the two types of breathing is illustrated in Figure 8. Two simulations were performed using typical ranges of pressure and volume changes for each type of breathing, see [West and Nadeau, 2003]. To do so, the reference stiffness K of the inner surface is tuned to 0.2 kPa. All the parameters used for each simulation are given in Table 4. Note that
under mechanical ventilation the pleural pressure has the nature of a reaction force, hence it is not prescribed. The total stress in the solid part $\Sigma_s$ has then two components: the first one is linked to the skeleton energy and describes the traction coming from the total loading, whereas the second one comes from the fluid pressure which leads to the solid part compression. These two components are presented separately in the figure.

## 3.5 Pulmonary symphysis

<table><caption>Table 5: Parameters used for the comparison between a normal case and a pathological case with pleural symphysis. Top rows correspond to the loading parameters, fluid pressure and pleural pressures, whereas bottom rows feature the coefficients of the constitutive behavior of the pulmonary mixture. The solid bulk modulus is not needed since $p_f = 0$.</caption><thead><tr><th>Name</th><th>Symbol</th><th>Unit</th><th>Healthy case</th><th>Symphysis case</th></tr></thead><tbody><tr><td>Air pressure</td><td>p<sub>f</sub></td><td>[kPa]</td><td>0</td><td>0</td></tr><tr><td>End-exhalation pleural pressure</td><td>p<sub>pl,e</sub></td><td>[kPa]</td><td>-0.5</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</td><td>-2</td></tr><tr><td>Sliding type</td><td></td><td></td><td>Frictionless</td><td>Tied</td></tr><tr><td rowspan="4">Skeleton energy parameters</td><td>&alpha;&#773;</td><td>[kPa]</td><td>8 &middot; 10<sup>-1</sup></td><td>8 &middot; 10<sup>-1</sup></td></tr><tr><td>&delta;</td><td>[-]</td><td>0.5</td><td>0.5</td></tr><tr><td>&beta;&#773;<sub>1</sub></td><td>[kPa]</td><td>0.1</td><td>0.1</td></tr><tr><td>&beta;&#773;<sub>2</sub></td><td>[kPa]</td><td>0.2</td><td>0.2</td></tr><tr><td>Solid bulk modulus</td><td>&kappa;&#773;</td><td>[kPa]</td><td>-</td><td>-</td></tr></tbody></table>

The proposed model also allows to study pathological lung conditions, for instance after a pleural symphysis. A pleural symphysis, where the two pleural membranes are joined together, can be performed to prevent pneumothorax from occurring or liquid from accumulating into the pleural space. Consequently, any sliding between the lungs and the thorax is prevented.

The modeling of a lung with pleural symphysis can be achieved using another type of contact – a tied contact – instead of the frictionless contact as previously used. A comparison between a normal case and this pathological case is now performed. The thorax is fixed, and the same amount of pressure as well as the same constitutive behavior are applied in both cases. The parameters used are given in Table 5.

The results are presented in Figure 9.

## 4 Discussion

We now discuss some specific points pertaining to the above results.

**Methods to enforce positive porosity** We considered two possible methods to enforce positive porosity, both in the direct and inverse problems. The first method, based on an additional pore energy was proposed by Chapelle and Moireau [2014], and is used here with a form of the pore
<figure><img src="image_11.png" /><figcaption>Figure 8: Solid stress components for both types of breathing, spontaneous and ventilated breathing. The quantities represented are the hydrostatic pressure of the solid Cauchy stress tensor, defined as $-\underline{\sigma}^H = \frac{1}{3} \text{tr } \underline{\sigma}$, associated with $\Sigma_{\underline{skel}}$ and $\Sigma_{\underline{p,s}}$. In the case of spontaneous breathing, the skeleton stress is in equilibrium with the mixture loading, here the pleural pressure of -0.8 kPa, and is negative reflecting the traction induced by the pleural pressure. In the case of ventilated breathing, the skeleton stress is smaller since the deformation is smaller, and $\underline{\sigma}_{\underline{p,s}}^H$ is positive, reflecting the compression induced by the fluid pressure imposed by the ventilator.</figcaption></figure>
energy that does not introduce residual stresses, and also enforces positive porosity in the inverse
poroelastic problem. The second method, based on a contact-like constraint, was introduced here.
In Figure 3 we showed the equivalence between the two methods in the limit where the pore energy
tends to zero.

**Positive porosity constraint** In the results presented in Fig. 6, we see that the estimation of the unloaded configuration induces a shift to the left of the porosity distribution – compared with the loaded configuration – because the porosity decreases when the lung deflates. When using no constraint, the left part of the distribution is negative, which is not physiological. However, the use of the contact-like formulation allows to keep the porosity positive everywhere. We can notice that the cells with vanishing porosity do not significantly impact the rest of the global distribution.

**Influence of the stiffness** We see in Fig. 7 how the stiffness parameter $\bar{\alpha}$ impacts the lung compliance. When $\bar{\alpha}$ increases, the volume change is lower for a given pressure, which means that the lung compliance decreases. Pathological cases could be represented in the model by adjusting the stiffness parameter. Indeed, the green curve ($\bar{\alpha} = 8 \cdot 10^{-1}$ kPa) could represent a fibrosis case, in which the lung becomes stiffer, the orange curve ($\bar{\alpha} = 8 \cdot 10^{-2}$ kPa) could stand for an healthy case, and the blue curve ($\bar{\alpha} = 8 \cdot 10^{-3}$ kPa) could be attributed to an emphysema case, where the lung is more compliant.

**Free vs. ventilated breathing** As the two components of the solid stress tensor do not represent the same physical effects, the total stress tensor $\Sigma_{\bar{s}}$ is not the relevant quantity to look at in order to study ventilation-induced injuries, which is why we represented the two components separately in Fig. 8. We see that the tensile stress component associated with skeleton deformation is somewhat

<figure><img src="image_12.png" /><figcaption>Figure 9: Comparison between two types of contact: the normal case with frictionless contact (left) and the case with pleural symphysis modeled with a tied contact (right). Both cases are shown in a coronal plane. The color scale represents the hydrostatic pressure of the solid Cauchy stress tensor (as defined in Section 3.3). The black line is the surface of the initial configuration.</figcaption></figure>
smaller in the ventilated case, indeed, whereas in this case we have a strong compressive component
due to internal fluid pressure, and therefore it it the latter component that is likely to be the cause
of VILI.

**Pleural symphysis** As seen in Fig. 9 the normal lung deforms mostly in the longitudinal direction since sliding is allowed against the thorax while the diaphragm pulls on the bottom part of the lung. By contrast, the lung with pleural symphysis cannot deform in the longitudinal direction to the same extent as the healthy lung because of symphysis. However, the pathological lung deforms more in the transverse direction. This behavior can explain the clinical observation that the rib cage diameter is larger for patients with pleural symphysis. Indeed, transverse deformation is then a compensation mechanism that allows to maintain an adequate amount of ventilation. Note that hydrostatic pressure in this case represents a scalar measure of the solid stress tensor – since porous pressure is zero – hence of the deformation tensor. We thus see in Fig. 9 that deformations are globally larger in the physiological case than in the pathological one, except near the internal surface, which was expected since boundary conditions associated with pulmonary symphysis are more constrained and the same loading is considered in both cases.

**Influence of gravity** Gravity has been included in our model equations for the sake of generality, but neglected in all our illustrations for practical purposes. It would require specific care, indeed, especially in the treatment of the boundary conditions associated with the unloaded configuration problem. Nevertheless, whereas taking gravity into account would make the model more physiological, it is expected that the resulting impact on pulmonary stresses and strains is of second order compared to breathing or ventilation pressure.

# 5 Conclusions

We have proposed a modeling framework based on large strain poromechanics to represent the behavior of lungs, and shown how the model can be personalized to individual cases when actual data – imaging, and possibly pleural pressure – are available. We have also illustrated the use of this model in a number of configurations motivated by clinical applications, including ventilated breathing.

In contrast to many pulmonary models focusing on air flows, this poromechanical model allows to study the parenchymal stress and strain with the use of a hyperelastic potential describing the skeleton behavior. The non-linear volume response of lungs to a pressure change as well as the solid quasi-incompressibility is well reproduced by the constitutive behavior.

Both effective and rescaled mechanical parameters can be used to parametrize the constitutive behavior, depending on whether or not porosity data are incorporated into the model. Unlike the effective parameters characterizing the mixture with the effect of the porosity, the rescaled parameters characterize the stiffness of the solid tissue. The effective parameters give the same type of porosity-dependent information as the clinical compliance, whereas the rescaled parameters bring an additional information which cannot be measured.

Since the constitutive behavior used in the model is nonlinear and the geometry comes from an *in vivo*, loaded configuration, the unloaded configuration needs to be retrieved. A special attention is paid to ensure positive porosity in the estimated unloaded configuration. We proposed two different
methods for this, including a novel contact-like material constraint that is the asymptotic limit of
the existing approach based on pore energy, and illustrated its response under various conditions.

We provided some preliminary illustrations of the effects of ventilated breathing, which actually
emphasizes the benefits of a poromechanical model compared with standard hyperelasticity, since
we can distinguish the very different effects induced in the solid stresses by the skeleton deformation
and by the fluid pressure. Based on this, some criteria could be proposed combining norms of both
quantities $\Sigma_{skel}$ and $\Sigma_{p,s}$ to quantify the risk of VILI. Moreover, as part of a poromechanical model,
stresses are macroscopic and homogenize the microscopic stresses, which are not available directly.
To access microscopic stress in the solid part, an alveolar geometry could be assumed and such
stresses could then be computed using the macroscopic deformation. Such microscopic information would provide the most relevant quantities for clinicians to study the impact of mechanical
ventilation on the pulmonary tissue.

One noteworthy limitation of the model is that it does not represent the dynamics of breathing, nor does it take into account the hysteresis behavior of the lung mechanical behavior that is commonly attributed to the pulmonary surfactant present on the alveoli surface [Wiechert, 2011]. However, these are not fundamental limitations, in the sense that the model could be extended to dynamics. Similarly, an energy describing the surfactant behavior could be incorporated into the solid free energy, in order to improve the model and make it more physiological.

In the future, this general pulmonary model should be validated against *in vivo* data and
could be applied in more complex and pathological cases involving mechanics – emphysema or
fibrosis for example – in which regional values of mechanical parameters could be defined instead
of homogeneous properties. Applied with a diseased patient, the personalized model could then
be used for diagnosis or classification purposes. Another question that such a model could help
investigate is the impact of drugs on the evolution of various diseases – by monitoring the model
parameters over the course of the disease evolution and treatments, in a patient specific approach
– in order to assist clinicians in deciding on the most appropriate medical treatment.

## Acknowledgements

We would like to thank Jean-Francois Bernaudin, Pierre-Yves Brillet, Hilario Nunes and Thomas Gille, clinicians from the Avicenne APHP Hospital, for their helpful discussions throughout the development of this work to understand lung physiology and pathologies. This work was partially supported by the French National Research Agency (ANR) under contracts numbers ANR-10-EQPX-37 and ANR-19-CE45-0007.

## References

Agostoni, E. (1972). Mechanics of the pleural space. *Physiological Reviews*, 52(1):57-128.

Al-Mayah, A., Moseley, J., Velec, M., and Brock, K. (2011). Toward efficient biomechanical-based deformable image registration of lungs for image-guided radiotherapy. *Physics in Medicine and Biology*, 56(15):4701-4713.

Alnæs, M., Blechta, J., Hake, J., Johansson, A., Kehlet, B., Logg, A., Richardson, C., Ring, J., Rognes, M. E., and Wells, G. N. (2015). The FEniCS Project Version 1.5. *Archive of Numerical Software*, 3(100).
Álvarez-Barrientos, F., Hurtado, D. E., and Genet, M. (2021). Pressure-driven micro-poromechanics: A variational framework for modeling the response of porous materials. *International Journal of Engineering Science*, 169:103586.

Baffico, L., Grandmont, C., and Maury, B. (2010). Multiscale Modeling of the Respiratory Tract. *Mathematical Models and Methods in Applied Sciences*, 20(01):59–93.

Baudet, V., Villard, P.-F., Jaillet, F., Beuve, M., and Shariat, B. (2003). Towards accurate tumour tracking in lungs. In *Proceedings on Seventh International Conference on Information Visualization, IV 2003.*, pages 338–343. IEEE Comput. Soc.

Berger, L., Bordas, R., Burrowes, K., Grau, V., Tavener, S., and Kay, D. (2016). A poroelastic model coupled to a fluid network with applications in lung modelling. *International Journal for Numerical Methods in Biomedical Engineering*, 32(1).

Biot, M. A. (1941). General Theory of Three-Dimensional Consolidation. *Journal of Applied Physics*, 12(2):155–164.

Biot, M. A. and Temple, G. (1972). Theory of Finite Deformations of Porous Solids. *Indiana University Mathematics Journal*, 21(7):597–620.

Bordas, R., Lefevre, C., Veeckmans, B., Pitt-Francis, J., Fetita, C., Brightling, C. E., Kay, D., Siddiqui, S., and Burrowes, K. S. (2015). Development and Analysis of Patient-Based Complete Conducting Airways Models. *Plos One*, 10(12).

Brock, K. K., Sharpe, M. B., Dawson, L. A., Kim, S. M., and Jaffray, D. A. (2005). Accuracy of finite element model-based multi-organ deformable image registration: Accuracy of FEM-based multi-organ deformable image registration. *Medical Physics*, 32(6Part1):1647–1659.

Budiansky, B. and Kimmel, E. (1987). Elastic Moduli of Lungs. *Journal of Applied Mechanics*, 54(2):351–358.

Burrowes, K. S., Iravani, A., and Kang, W. (2019). Integrated lung tissue mechanics one piece at a time: Computational modeling across the scales of biology. *Clinical Biomechanics*, 66:20–31.

Burtschell, B. (2016). *Mechanical Modeling and Numerical Methods for Poromechanics: Application to Myocardium Perfusion*. PhD thesis, École Polytechnique, Palaiseau, France.

Caruel, M., Chabiniok, R., Moireau, P., Lecarpentier, Y., and Chapelle, D. (2014). Dimensional reductions of a cardiac model for effective validation and calibration. *Biomechanics and Modeling in Mechanobiology*, 13(4):897–914.

Cazeaux, P. and Grandmont, C. (2015). Homogenization of a multiscale viscoelastic model with nonlocal damping, application to the human lungs. *Mathematical Models and Methods in Applied Sciences*, 25(06):1125–1177.

Chabiniok, R., Wang, V. Y., Hadjicharalambous, M., Asner, L., Lee, J., Sermesant, M., Kuhl, E., Young, A. A., Moireau, P., Nash, M. P., Chapelle, D., and Nordsletten, D. A. (2016). Multiphysics and multiscale modelling, data–model fusion and integration of organ physiology in the clinic: Ventricular cardiac mechanics. *Interface Focus*, 6(2):20150083.
Chapelle, D., Gerbeau, J.-F., Sainte-Marie, J., and Vignon-Clementel, I. E. (2010). A poroelastic model valid in large strains with applications to perfusion in cardiac modeling. *Computational Mechanics*, 46(1):91–101.

Chapelle, D. and Moireau, P. (2014). General coupling of porous flows and hyperelastic formulations—From thermodynamics principles to energy balance and compatible time schemes. *European Journal of Mechanics - B/Fluids*, 46:82–96.

Ciarlet, P. and Geymonat, G. (1982). Sur les lois de comportement en élasticité non-linéaire compressible. *C. R. Acad. Sci. Paris Sér. II*, 295:423–426.

Clark, A. R., Kumar, H., and Burrowes, K. (2017). Capturing complexity in pulmonary system modelling. *Proceedings of the Institution of Mechanical Engineers, Part H: Journal of Engineering in Medicine*, 231(5):355–368.

Coussy, O. (2004). *Poromechanics*. Wiley, Chichester.

Dale, P., Matthews, F. L., and Schroter, R. C. (1980). Finite element analysis of lung alveolus. *Journal of Biomechanics*, 13(10):865–873.

Demiray, H. (1972). A note on the elasticity of soft biological tissues. *Journal of Biomechanics*, 5(3):309–311.

Denny, E. and Schroter, R. (2006). A model of non-uniform lung parenchyma distortion. *Journal of Biomechanics*, 39(4):652–663.

Dormieux, L., Kondo, D., and Ulm, F.-J. (2006). *Microporomechanics*. John Wiley & Sons, Chichester, West Sussex, England ; Hoboken, NJ.

Eom, J., Xu, X. G., De, S., and Shi, C. (2010). Predictive modeling of lung motion over the entire respiratory cycle using measured pressure-volume data, 4DCT images, and finite-element analysis: A novel physics-based respiratory modeling for entire breathing cycle. *Medical Physics*, 37(8):4389–4400.

Fuerst, B., Mansi, T., Carnis, F., Salzle, M., Zhang, J., Declerck, J., Boettger, T., Bayouth, J., Navab, N., and Kamen, A. (2015). Patient-Specific Biomechanical Model for the Prediction of Lung Motion From 4-D CT Images. *IEEE Transactions on Medical Imaging*, 34(2):599–607.

Fung, Y. C. (1981). *Biomechanics*. Springer New York, New York, NY.

Gee, M. W., Förster, C., and Wall, W. A. (2010). A computational strategy for prestressing patient-specific biomechanical problems under finite deformation. *International Journal for Numerical Methods in Biomedical Engineering*, 26(1):52–72.

Genet, M. (2019). A relaxed growth modeling framework for controlling growth-induced residual stresses. *Clinical Biomechanics*, 70:270–277.

Genet, M., Chuan Lee, L., Ge, L., Acevedo-Bolton, G., Jeung, N., Martin, A., Cambronero, N., Boyle, A., Yeghiazarians, Y., Kozerke, S., and Guccione, J. M. (2015a). A Novel Method for Quantifying Smooth Regional Variations in Myocardial Contractility Within an Infarcted Human Left Ventricle Based on Delay-Enhanced Magnetic Resonance Imaging. *Journal of Biomechanical Engineering*, 137(8):081009.
Genet, M., Rausch, M., Lee, L., Choy, S., Zhao, X., Kassab, G., Kozerke, S., Guccione, J., and Kuhl, E. (2015b). Heterogeneous growth-induced prestrain in the heart. *Journal of Biomechanics*, 48(10):2080-2089.

Gibson, G. and Pride, N. (1976). Lung distensibility. The static pressure-volume curve of the lungs and its use in clinical assessment. *British Journal of Diseases of the Chest*, 70:143-184.

Goligher, E. C., Ferguson, N. D., and Brochard, L. J. (2016). Clinical challenges in mechanical ventilation. *Lancet (London, England)*, 387(10030):1856-1866.

Govindjee, S. and Mihalic, P. A. (1998). Computational methods for inverse deformations in quasi-incompressible finite elasticity. *International Journal for Numerical Methods in Engineering*, 43(5):821-838.

Han, L., Dong, H., McClelland, J. R., Han, L., Hawkes, D. J., and Barratt, D. C. (2017). A hybrid patient-specific biomechanical model based image registration method for the motion estimation of lungs. *Medical Image Analysis*, 39:87-100.

Hinz, B. and Suki, B. (2016). Does Breathing Amplify Fibrosis? *American Journal of Respiratory and Critical Care Medicine*, 194(1):9-11.

Hurtado, D. E., Villarroel, N., Andrade, C., Retamal, J., Bugedo, G., and Bruhn, A. (2017). Spatial patterns and frequency distributions of regional deformation in the healthy human lung. *Biomechanics and Modeling in Mechanobiology*, 16(4):1413-1423.

Im, J. G., Webb, W. R., Rosen, A., and Gamsu, G. (1989). Costal pleura: Appearances at high-resolution CT. *Radiology*, 171(1):125-131.

Kallet, R. H. (2015). A Comprehensive Review of Prone Position in ARDS. *Respiratory Care*, 60(11):1660-1687.

Kaul, H. (2019). Respiratory healthcare by design: Computational approaches bringing respiratory precision and personalised medicine closer to bedside. *Morphologie*, 103(343):194-202.

Kowalczyk, P. (1993). Mechanical model of lung parenchyma as a two-phase porous medium. *Transport in Porous Media*, 11(3):281-295.

Lai-Fook, S. J. (2004). Pleural Mechanics and Fluid Exchange. *Physiological Reviews*, 84(2):385-410.

Le Tallec, P. (1994). Numerical methods for nonlinear three-dimensional elasticity. In Ciarlet, P. G. and Lions, J.-L., editors, *Handbook of Numerical Analysis*, volume 3. Elsevier.

Lee, L. C., Genet, M., Dang, A. B., Ge, L., Guccione, J. M., and Ratcliffe, M. B. (2014). Applications of computational modeling in cardiac surgery. *Journal of Cardiac Surgery*, 29(3):293-302.

Logg, A., Mardal, K., and Wells, G. (2012). *Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book*. Lecture Notes in Computational Science and Engineering. Springer Berlin Heidelberg.

Millar, A. B. and Denison, D. M. (1989). Vertical gradients of lung density in healthy supine men. *Thorax*, 44(6):485-490.
Miserocchi, G. (1997). Physiology and pathophysiology of pleural fluid turnover. *European Respiratory Journal*, 10(1):219-225.

Noppen, M., De Waele, M., Li, R., Gucht, K. V., D'Haese, J., Gerlo, E., and Vincken, W. (2000). Volume and Cellular Content of Normal Pleural Fluid in Humans Examined by Pleural Lavage. *American Journal of Respiratory and Critical Care Medicine*, 162(3):1023-1026.

Ogden, R. W. (1972). Large deformation isotropic elasticity: On the correlation of theory and experiment for compressible rubberlike solids. *Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences*, 328(1575):567-583.

OMS (2016). The top 10 causes of death in 2016. http://www.who.int/news-room/fact-sheets/detail/the-top-10-causes-of-death.

Patte, C., Brillet, P.-Y., Fetita, C., Gille, T., Bernaudin, J.-F., Nunes, H., Chapelle, D., and Genet, M. (2020). Estimation of patient-specific mechanical parameters in pulmonary diseases. *VPH2020 conference*.

Rausch, M. K., Genet, M., and Humphrey, J. D. (2017). An augmented iterative method for identifying a stress-free reference configuration in image-based biomechanical modeling. *Journal of Biomechanics*, 58:227-231.

Reeve, A. M., Nash, M. P., Taberner, A. J., and Nielsen, P. M. F. (2014). Constitutive Relations for Pressure-Driven Stiffening in Poroelastic Tissues. *Journal of Biomechanical Engineering*, 136(8):081011.

Richardson, S., Gamage, T. P. B., HajiRassouliha, A., Jackson, T., Hedges, K., Clark, A., Taberner, A., Tawhai, M. H., and Nielsen, P. M. F. (2019). Towards a Real-Time Full-Field Stereoscopic Imaging System for Tracking Lung Surface Deformation Under Pressure Controlled Ventilation. In Nielsen, P. M. F., Wittek, A., Miller, K., Doyle, B., Joldes, G. R., and Nash, M. P., editors, *Computational Biomechanics for Medicine*, pages 119-130. Springer International Publishing, Cham.

Rivlin, R. S. (1948). Large elastic deformations of isotropic materials. I. Fundamental concepts. *Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences*, 240(822):459-490.

Sarabia-Vallejos, M. A., Zuñiga, M., and Hurtado, D. E. (2019). The role of three-dimensionality and alveolar pressure in the distribution and amplification of alveolar stresses. *Scientific Reports*, 9(1):8783.

Scheiner, S., Pivonka, P., and Hellmich, C. (2013). Poromechanical Stimulation of Bone Remodeling: A Continuum Micromechanics-Based Mathematical Model and Experimental Validation. In *Poromechanics V*, pages 1867-1876, Vienna, Austria. American Society of Civil Engineers.

Sellier, M. (2011). An iterative method for the inverse elasto-static problem. *Journal of Fluids and Structures*, 27(8):1461-1470.

Seyfi, B., Santhanam, A. P., and Ilegbusi, O. J. (2016). A Biomechanical Model of Human Lung Deformation Utilizing Patient-Specific Elastic Property. *Journal of Cancer Therapy*, 07(06):402-415.
Sundaram, T. A. and Gee, J. C. (2005). Towards a model of lung biomechanics: Pulmonary kinematics via registration of serial lung images. *Medical Image Analysis*, 9(6):524-537.

Tawhai, M. H., Nash, M. P., Lin, C.-L., and Hoffman, E. A. (2009). Supine and prone differences in regional lung density and pleural pressure gradients in the human lung with constant shape. *Journal of Applied Physiology*, 107(3):912-920.

Wang, N.-S. (1975). The Preformed Stomas Connecting the Pleural Cavity and the Lymphatics in the Parietal Pleura. *American Review of Respiratory Disease*, 111(1).

Weiss, J. A., Maker, B. N., and Govindjee, S. (1996). Finite element implementation of incompressible, transversely isotropic hyperelasticity. *Computer Methods in Applied Mechanics and Engineering*, 135(1-2):107-128.

Werner, R., Ehrhardt, J., Schmidt, R., and Handels, H. (2009). Patient-specific finite element modeling of respiratory lung motion using 4D CT image data: Finite element modeling of respiratory lung motion. *Medical Physics*, 36(5):1500-1511.

West, J. B. and Matthews, F. L. (1972). Stresses, strains, and surface pressures in the lung caused by its weight. *Journal of Applied Physiology*, 32(3):332-345.

West, J. B. and Nadeau, M. (2003). *La physiologie respiratoire*.

Wiechert, L. (2011). *Computational Modeling of Multi-Field and Multi-Scale Phenomena in Respiratory Mechanics*. PhD thesis, TU Munich, Germany.

Wiechert, L. and Wall, W. (2010). A nested dynamic multi-scale approach for 3D problems accounting for micro-scale multi-physics. *Computer Methods in Applied Mechanics and Engineering*, 199(21-22):1342-1351.

Winslow, R. L., Trayanova, N., Geman, D., and Miller, M. I. (2012). Computational Medicine: Translating Models to Clinical Care. *Science Translational Medicine*, 4(158):158rv11.

Zeng, Y. J., Yager, D., and Fung, Y. C. (1987). Measurement of the Mechanical Properties of the Human Lung Tissue. *Journal of Biomechanical Engineering*, 109(2):169-174.

Zhang, T., Orton, N. P., Mackie, T. R., and Paliwal, B. R. (2004). Technical note: A novel boundary condition using contact elements for finite element based deformable image registration. *Medical Physics*, 31(9):2412-2415.

Zocchi, L. (2002). Physiology and pathophysiology of pleural fluid turnover. *European Respiratory Journal*, 20(6):1545-1558.