Pressure- and time-dependent alveolar recruitment/derecruitment in a spatially resolved patient-specific computational model for injured human lungs

We present a novel computational model for the dynamics of alveolar recruitment/derecruitment (RD), which reproduces the underlying characteristics typically observed in injured lungs. The basic idea is a pressure- and time-dependent variation of the stress-free reference volume in reduced dimensional viscoelastic elements representing the acinar tissue. We choose a variable reference volume triggered by critical opening and closing pressures in a time-dependent manner from a straightforward mechanical point of view. In the case of (partially and progressively) collapsing alveolar structures, the volume available for expansion during breathing reduces and vice versa, eventually enabling consideration of alveolar collapse and reopening in our model. We further introduce a method for patient-specific determination of the underlying critical parameters of the new alveolar RD dynamics when integrated into the tissue elements, referred to as terminal units, of a spatially resolved physics-based lung model that simulates the human respiratory system in an anatomically correct manner. Relevant patient-specific parameters of the terminal units are herein determined based on medical image data and the macromechanical behavior of the lung during artificial ventilation. We test the whole modeling approach for a real-life scenario by applying it to the clinical data of a mechanically ventilated patient. The generated lung model is capable of reproducing clinical measurements such as tidal volume and pleural pressure during various ventilation maneuvers. We conclude that this new model is an important step toward personalized treatment of ARDS patients by considering potentially harmful mechanisms - such as cyclic RD and overdistension - and might help in the development of relevant protective ventilation strategies to reduce ventilator-induced lung injury (VILI).


INTRODUCTION
In patients with acute respiratory distress syndrome (ARDS), mechanical ventilation is a potentially life-saving treatment.Nevertheless, improperly adjusted ventilator settings can lead to ventilator-associated lung injury (VILI) 1,2,3 .Despite intense research in this field, which has among other things revealed two major contributors to VILI -i.e., overdistension (volutrauma) and cyclic (re-)opening (atelectrauma) of lung structures 3,4 -, the development of protective ventilation strategies stagnated over the last two decades 5,6 .
The main reason for this halt is the inaccessibility of insights into regional mechanics of patients' lungs and damaging processes occurring during ventilation 7 .The situation is particularly exacerbated by the unique heterogeneous pathology of every diseased organ causing irregular and unpredictable tissue straining, tidal recruitment and distribution of ventilation 7,8,9,10 .Local stress raisers, that is, sites of high and thus harmful stress, which presumably occur in the presence of lung inhomogeneity and dysfunction 11,12 and were identified as the origin of injuries 4,10,13 , cannot be determined in routine clinical practice, i.e., the impact of an applied ventilation protocol on the regional microscale -may it be beneficial or harmful -is hardly assessable, in particular for a specific patient 6,7,8,9 .As a result, the current clinical practice of generic ventilation management and ventilator adjustment at the bedside 14 reaches its limits.An individual heterogeneous lung injury needs an individual therapy management 15,16 .However, such management is only possible by achieving a locally deeper, patient-specific understanding of ARDS lung mechanics 7,15 .
Computational models of the human lung are promising tools for serving this purpose non-invasively.A broad range of modelling approaches already exist 17,18 .By virtue of their high computational efficiency, reduced dimensional models in particular are feasible for use in clinical application.
Many models that explicitly address both RD and OD are of a single-compartmental and phenomenological nature 27,29,31 .They allow for conclusions about general correlations, e.g., the categorization of lung injury by mechanical power dissipation and strain heterogeneity 29 , or, that volu-and atelectrauma contribute to VILI in a combined manner and not independently 27 .Since neither local effects nor the patient-specific pathology are captured by these models, but these factors gain in importance when needing to adapt mechanical ventilation adequately 7,32,33 , their benefit to clinical therapy management is limited.
In addition, the few multi-compartment models of the lung mimicking RD and OD have shortcomings, e.g., due to not considering local pathology and neglecting time dependence of RD 19,21,28 .Especially the latter is known as a relevant phenomenon of RD 34,35 , and recent studies have demonstrated that timing has a great impact on the ventilation of injured lungs 4,6,10,32,36 .In a previous work 37 , we included the individual regional heterogeneity into a physics-based, anatomically accurate reduced dimensional lung model by incorporating time-dependent RD dynamics 34,38 on the conducting airway level and linking them to anatomical and pathological specifications.It was a first step toward including the potential for atelectrauma (however without explicity addressing and investigating it) as well as volutrauma by locating and characterizing mechanical stress foci more accurately.Nevertheless, the research revealed the challenge in determining RD parameters uniquely for a patient, considering that physical conditions are not accessible locally and might differ regionally in the organ, and also from one patient to another.The model also neglected the probable presence of alveolar RD, which can result in an inaccurate regional volume capacity and, thus, deficient estimation of local straining on the one hand, and the cause of atelectrauma in the lung parenchyma on the other hand.
The current objective in this paper is to overcome these issues and advance the personalization of computational lung models by providing a tool that can be used to assess the VILI potential of certain ventilation profiles for a patient.We introduce a novel approach to modeling pressure-and time-dependent alveolar RD dynamics integrated into a viscoelastic component representing pulmonary tissue.One single (de-)recruitable viscoelastic component can represent both alveolar distension and RD.The idea being presented is motivated by previous studies 27,29,31,35 , but it is based more on mechanical than phenomenological principles via the concept of modifying the stress-free reference volume of tissue in order to mimic RD, and it is meant to be applied as a multi-compartment model.This new model for RD is employed in the tissue elements (terminal units) of our comprehensive lung model.To enable patient-specific model calibration, we present a generic method that tailors the numerous model parameters (especially those related to RD) to ARDS patients based on ventilation characteristics and on their pathology extracted from medical image data.By applying this procedure to a specific patient, we examine the model's capacity to reproduce the clinical ventilation protocol performed at the bedside, and, thus, its predictive capability for this patient.Further, we investigate local straining and RD, the two mechanisms predominantly responsible for VILI.In the long run, such a digital twin can help to evaluate the injury potential of ventilation protocols for each patient specifically, and to eventually minimize VILI individually.

Physiological background
Lung tissue is a very delicate structure when viewing the small airways, alveolar ducts, and alveoli.The actual state of diseased tissue in vivo is very complex and remains difficult to resolve at the micro level.It can exhibit air-less and potentially recruitable lung units, air spaces flooded with inflammatory fluid, abnormal swelling of the alveolar wall, or a combination of these 39 .Medical imaging enables identification of these scattered pathological regions, but their exact condition cannot be determined 39 .
Derecruitment due to lung edema may reduce the capacity of the tissue to distend because of fluid occupation or stiffening of alveolar walls, and can be regarded accordingly when modeling a patient's lung (e.g., 26 ).In contrast to the relatively constant condition of edema presence, the phenomenon of cyclic intra-tidal opening and closing of lung structures has a very dynamic nature, which has been researched intensely.Ghadiali and Huang 40 presented a comprehensive review on findings about RD and its impact on VILI.We will therefore only recapitulate herein the main characteristics of alveolar RD found in the literature and forming the basis of our modeling approach and underlying assumptions.
Alveolar ducts, and therefore alveolar tissue, change in size not only due to straining, but also because of (de-)folding and closure or opening 41 .The collapse of injured distal airspaces happens heterogeneously 41 , so not all at once on the acinar level.This (gradual) derecruitment due to closure of lung units reduces the amount of tissue available for distension during a breath cycle, thus resulting -precisely as in edema formation -in a decrease in lung compliance (meaning the change in lung volume per change in transpulmonary pressure) 39,42,43 .Conversely, recruitment brings about an increase in lung compliance.
To the extent evident, RD is primarily triggered by pressure, but exhibits a pronounced time-dependent behavior 4,10,34,35,44 .The pressure and time dependence each have individual characteristics: • The prevailing pressure in the lung defines the eventual amount of recruited volume in a lung (in the subsequent paragraph denoted as  f inal ) 35 .Further, the opening pressure of a lung unit is typically higher than the pressure at which a fully open and thus stable alveolar structure is supposed to close 9,19,27,38,45 .Viewing RD in an injured lung, the critical opening and closing pressures are not necessarily statistically distributed, but significantly scattered 10 .
• The time dependence of RD has been observed in various experimental and clinical setups 4,10,35,46,47,48 .The volume increase over time due to recruitment has often been described by a relaxation relationship following 46,47,48 , with the final volume  f inal approximated over time during constant ventilation pressure and the time constant  specifying the speed of opening toward  f inal .Previous studies 46,47,48 employed this equation to specify the chronological progression of globally observed RD in ARDS lungs.Using the same relation in our model (see Section 2.2) provides a valuable indication for a reasonable choice of  when tailoring the model parameters to a specific patient.Further, Albert et al. 35 showed in mice with saline lavaged lungs that the major alveolar recruitment happens after the first two seconds, regardless of the prevailing pressure level, i.e., a difference in the recruitment pressure level did not influence the time constant .
Remark: In a diseased lung, there are supposedly two different types of time dependencies, i.e., viscous effects that originate from normal tissue expansion and those that arise due to RD 29,48 .Our model accounts for both types of time dependencies.

Modeling of alveolar recruitment/derecruitment
We use a well-established model element for acinar tissue 26,49,50 and enhance it with a novel, mechanically motivated model for alveolar RD based on the RD characteristics presented, which we will describe in greater detail hereinafter.Given that this model is eventually used at the outlets of conducting reduced dimensional airway elements (see underlying model in Section A in the Appendix) and mimics the subsequent acinar region, we will hereinafter refer to this tissue element as the terminal unit.The viscoelastic behavior of lung parenchyma can be described in a reduced-dimensional manner by the arrangement of springs and dashpots (generalized Maxwell model 51 ).In this study, we assume every terminal unit to consist of  ad alveolar ducts (Figure 1), each represented by a Kelvin-Voigt body (spring and dashpot connected in parallel) in parallel with a Maxwell body (spring and dashpot connected in series) 26,37,49,50 .A terminal unit then follows where  is the gas flow into the terminal unit,  is the pressure difference between the alveolar pressure  alv inside the terminal unit and the surrounding pressure  pl (see Section 3.3 for more details),  and  a are linear dashpots modeling time-dependent effects of the viscoelastic tissue distension, and  E1 and  2 are the non-linear and linear springs, respectively.Especially in the context of OD, mimicking the non-linear behavior of lung tissue is a crucial modeling aspect.This effect is achieved by the non-linear spring  E1 which represents the static pressure-volume relationship of a terminal unit due to its arrangement in the generalized Maxwell model.The non-linear pressure-volume behavior of  E1 is derived from a purely volumetric deformation of an Ogden-type material 26,37,52 yielding where  denotes the current gas volume of a terminal unit and  0 the reference value of  in the stress-free state. and  are slope-and curvature-shaping parameters, respectively.In the case of (partially and progressively) collapsing alveolar structures or the infiltration of inflammatory liquid and edema, we assume that the alveolar volume available for expansion reduces, as outlined above.Similarly, the recruitment of air spaces increases the available volume.These changes in volume eventually decrease or increase the overall compliance of a terminal unit, respectively.Considering both this effect and the time and pressure dependency of RD, we included RD dynamics in the model of a terminal unit presented above by making  0 in Eq. (2) dependent on pressure and time (Figure 2), which is referred to as the current reference volume  0 ( , ).

Pressure dependence
The  experienced by a terminal unit regulates  0 as follows: Below a minimal critical pressure  crit,min , the reference gas volume tends to fully collapse and reach its minimal volume  0,min (Figure 2, left).Above  crit,min , the target value of the reference gas volume linearly depends on  until it reaches the maximal critical pressure  crit,max =  crit,min + Δ max−min .Above  crit,max ,  0 approaches the maximal reference gas volume  0,max , indicating a fully open state of all alveoli.The relation between  0 and its V 0,targ (P) with where the pressure dependent reference volume calculated from Eq. ( 3) is denoted as target reference volume  0,targ ( ), to which  0 ( , ) tends in a time dependent manner.The time dependence of  0 ( , ) is described in the following paragraph.This approach makes it possible to model not only an open and closed state in a binary manner 37 when using a RD model in the conducting airways, but also the gradual RD and partial collapse of the heterogeneously (de-)recruiting alveolar structure represented by one terminal unit.In this state, the open tissue proportion can still expand and potentially also experience OD.The current reference volume  0 ( , ) also provides an indication of the opening degree of a terminal unit.To supplement the RD model by also considering that critical opening pressures usually exceed the closing pressures of a stable or recruited alveolar structure (see Section 2.1), we introduce two different paths for the recruitment and derecruitment process of a terminal unit, which differ by a shift Δ op−cl of the corresponding critical pressures (see opening and closing path in Figure 2, left).
Due to the lack of knowledge on the shape of the pressure-volume relationship for opening of a conglomerate of alveoli, the assumption of a linear relationship between reference gas volume and applied pressure  is a first starting point.To the best of our knowledge, the literature on the underlying characteristics is very sparse.Thus, this approach can undoubtedly be refined in future work as more insights become available from experiments or resolved simulations.

Time dependence
As introduced in Section 2.1, we used a time dependence relationship previously applied in a variety of contexts in order to describe the temporal alveolar (de-)recruitment in injured lungs:  0 ( , ) yields  0,targ ( ) over time according to where  0,init is an initial reference gas volume and  is a time constant specifying the delay in RD of the terminal unit.Since  0 ( , ) tends to  0,targ ( ) without ever reaching it mathematically, we introduce the tolerance factor  V 0 to ensure that a terminal unit can be opened or closed entirely.If a previously closed terminal unit (i.e., with a  0,targ ( ) moving on the opening path, see Figure 2, left) exceeds ( 1 −  V 0 )  0,max , it is assumed to fully recruit, and the opposite is assumed when falling below )  0,min .Such a change of state eventually triggers the switch between the opening and closing paths and, therefore, between the relevant critical pressures.
Due to their clustered structure, we assume that the  ad alveolar ducts of a single terminal unit in Eq. ( 1) have similar RD parameters used in Eqs. ( 3) -( 5) modeled by a linear relation between  and  0,targ for the entire terminal unit.This approach is based on the assumption that the subregion of the respiratory zone modeled by a terminal unit is governed by a specific, relatively uniform (patho-) physiological condition.One method used to determine the parameters for a single terminal unit based on patient imaging and ventilation data is described in Section 3.4, resulting in individual parameter sets for each element in the framework of a multi-compartment patient-specific lung model.
As described in Section 2.1, alveolar derecruitment might also occur due to infiltration or edema in the alveolar tissue.This occupation of air spaces eventually decreases the distensible volume.In this study, we account for the accumulated fluid in the RD model by a constant amount of volume by which  0,max is reduced (Section 3.4).
We implemented the presented model by numerically discretizing Eqs. ( 1) -( 5) with the first-order Euler scheme.Similar to Bates and Irvin 34 , we performed several tests in an idealized setup to confirm that the novel RD model for terminal units can reproduce typical RD characteristics, and evaluated the influence of individual model parameters of Eqs.(3) - (5).However, in order to not overload this paper we refrain from showing the underlying tests.

SETTING UP A PATIENT-SPECIFIC MODEL
In the long run, one objective of the proposed alveolar RD model is that of being applicable to patient-specific comprehensive multi-compartment lung models with spatial resolution.In addition to the anatomically accurate geometry generation, this objective also requires the ability to realistically tailor all underlying parameters to a patient, especially the numerous RD parameters for each terminal unit.The parametrization should be feasible with a reasonable effort and include the ability to capture regional heterogeneity as observable in medical image data.For this purpose, we developed a generic procedure that considers a patient's pathology, as extracted from medical imaging, and the mechanical behavior of the patient's respiratory system obtained from a few clinical ventilation measurements.In the following, we describe the full approach to generating and calibrating an anatomically accurate, comprehensive multi-compartment lung model for a specific patient incorporating the presented model for alveolar RD.

Clinical data
Image data Generating the model of a patient-specific lung requires a three-dimensional thoracic CT scan of the ventilated patient recorded at a known level of -usually positive -end-expiratory pressure (PEEP) along normal ventilation.Such a CT scan is usually part of the standard protocol for ARDS patients when admitted to the ICU.Based on this scan, we extract the geometry and individual pathology of the lung indicated by gray values of the image voxels.Ventilation data To calibrate the model parameters, we use bedside ventilation measurements, including the pressure at the airway opening, the tracheal airflow entering the lung, and the esophageal pressure as a surrogate for the pleural pressure.We specifically also use measurement of a normal breath cycle of the ventilation mode during which the CT scan has been recorded, and a quasi-static inflation maneuver to obtain information about the mechanical properties of the patient's respiratory system.Transpulmonary pressures used in the following are computed as the difference between the pressure at the airway opening and the pleural pressure.

Geometry generation
The conducting airway tree of a lung model starts from the centerline of the tracheal tree which reaches to the lobar bronchi leading into the individual lung lobes.The geometry and dimension of all these components are extracted from the CT scan.Due to the limited resolution of the medical images, we apply a space-filling tree growing algorithm for the regions beyond the lobar bronchi 53 , and extended and used in several of our previous studies 26,49,50 .A single airway segment of the conducting airway tree is modelled as described in Section A in the Appendix.Every terminal airway of the tree supplies gas to a terminal unit attached to it (see underlying model in Section 2.2).
To assign a realistic volume and Hounsfield Unit (HU) to each terminal unit, we slightly extended the tree growing procedure: The cloud of voxels constituting the lung region in the CT scan is continuously split into subgroups with proceeding branching of the airway tree.All voxels which remain at a terminal airway are assigned to the adherent terminal unit (Figure 3).We can in this way attribute a total volume and a mean HU -calculated from the voxels -to each terminal unit.This information enables characterization of the composition of the total volume of gas ( PEEP ) and water ( water ) following HU = -1000 equals air and HU = 0 equals water (∼ tissue and edema liquid) 39 .These quantities are used to parametrize the novel RD model in Section 3.4.

FIGURE 3
Conducting airway tree of a computational lung model resulting from the space-filling tree growing algorithm, illustrated using a slice from the underlying CT scan (left); the voxel cloud remaining at a terminal airway after continuous splitting during the tree growing algorithm (right) defines the water and gas volume and the mean HU of a terminal unit.

Individualized pressure boundary conditions
Building on our previous lung model 37 , but with slight detail modifications, the external pressure  pl acting on the terminal units in simulations consists of two components: the variable, volume-dependent pressure  vol pl and the static contribution  weight pl introducing the weight of the lung regions above the units in question as an additional super-imposed pressure, and thus depending on the height (Figure 4) 54,55 .Hence, being the sum of these two constituents, the pleural pressure reads  pl =  vol pl +  weight pl .Basically,  vol pl is the passive pressure exerted by the sedated chest wall of a patient against the expanding lung.In contrast to our previous work 37 , we use a non-linear progression for the relationship between lung volume and pleural pressure reading with where  v ,  v ,  v and  v are determined for the individual patient. f rac is the volume share of the volume increase in all terminal units,  tot −  tot,PEEP , from the overall air volume in the lung at PEEP,  tot,PEEP , calculated from the CT scan, and the increase of air volume at end-inspiration of the measured quasi-static pressure-volume curve of the patient,  tot,max −  tot,PEEP .
In diseased lungs, the pressure gradient occurring across the lung due to its own weight is often more pronounced due to additional weight from pathological fluid accumulations in the organ 54 .For determination of the patient-specific  weight pl , which is in the patient's supine position a function of the ventral-to-dorsal height of the lung, ℎ, we largely followed the approach used by Pelosi et al. 54 .Instead of dividing the lung only into ten vertical intervals, however, we discretize the lung in voxel slices to calculate the respective superimposed pressure in each layer.For the sake of precision, we fit the resulting discrete pressure values to a quartic relationship in place of a quadratic curve 54 , yielding where  w ,  w ,  w ,  w and  w are fitting parameters.To enforce  weight pl = 0 at the reference point of measurement of  pl , which was made with the esophageal balloon in the esophagus (see positive and negative range of  weight pl in Figure 4), we amended Eq. ( 8) with the height of the measurement spot, ℎ balloon , determined from the CT scan 56 .

Image-and ventilation-based parametrization of terminal units
In order to apply the presented model of alveolar RD to a real ARDS lung and individualize it for a patient, we propose a sophisticated approach to calibrate the numerous model parameters.It bases on a deep understanding of the model and lung mechanics, and uses image data and deliberately selected ventilation maneuvers to extract certain patient-specific information.In this way, we manage to handle the underlying complexity when applying the model to a real case and establish an innovative way to determine the model parameters for the respiratory mechanics and lung pathology of a specific patient.For the sake of completeness, before going into more detail about the calibration method, we would like to point out that there are less refined and not explicitly designed for our purpose, but still effective procedures for model parametrization, e.g., Bayesian inverse analysis methods, which our group is currently developing also in the context of biomechanical problems 57,58,59 .These methods might be of particular interest when using data from further or different measurement sources that are not used in the calibration method described below.
In the following, we propose a two layered procedure that parametrizes the terminal units: On the one hand, we determine the RD parameters in Eqs. ( 3) -( 5) for each terminal unit by means of an algorithm (Section 3.4.1)based on the mean gray value and the gas volume of the terminal unit assigned in the geometry generation, and based on its approximated height-dependent pressure state at end-inspiration and end-expiration during a measured normal ventilation breath cycle.The gray values herein indicate the pathological state of the terminal unit and thus enable adequate selection of its RD parameters to reproduce the regional heterogeneity of the lung.On the other hand, the aforementioned algorithm requires a set of input parameters (see a list of the input parameters in Table B2), which are optimized to match the overall pressure-volume behavior of all terminal units to the measurements of a normal breath cycle and a quasi-static inflation maneuver (Section 3.4.2).

Assumptions
To simplify the algorithm for RD parametrization and the optimization of its necessary input parameters, both of which are described in greater detail hereinafter, we for now base it on the following assumptions below.It is important to note that these assumptions are made only during model parametrization and do not hold for the whole computational model.
• We disregard the effects of the upstream airway tree, e.g., the pressure drop across the conducting airways due to their flow resistance.We thus assume that the pressure at the airway opening directly acts upon the terminal units, especially at quasi-static points during ventilation like end-inspiration, end-expiration and during quasi-static inflation, where this assumption approximately holds true.
• We only consider the elastic component of each terminal unit, i.e., the non-linear spring  E1 following the relation in Eq. ( 2) that incorporates RD dynamics, but disregard any effects of the spring  2 and the dashpots  and  a in Eq. ( 1).
• The pressure state within the lung during recording of the CT scan,  tp,PEEP , is PEEP −  pl , with the latter depending from  tot,PEEP and ℎ as depicted in Figure 4 and described in Section 3.3.In other words, a specific and identifiable heightdependent  tp,PEEP is assumed for each terminal unit at PEEP in the moment of capturing the CT scan which is relevant when determining the critical RD pressures 55,60 .
• We distinguish among three types of terminal units, depending on the HU assigned, as is similar to the assumed pathology condition in an ARDS lung 6,39,42 : normally aerated (including hyper-ventilated) terminal units with HU ≤ −500, poorly ventilated terminal units with −500 < HU < −100, and non-aerated terminal units with HU ≥ −100.
• We assume  =  insp =  exp of a terminal unit for the sake of simplicity in contrast to the model proposed by Bates and Irvin 34 in which airways tended to close faster than they opened.
• For each collapsed and poorly ventilated terminal unit, we assume a predefined constant ratio  edema of edema volume and the thereby filled tissue volume.Based on the volumes  water and  PEEP extracted from the CT scan, we assume for  tissue that it is comprised of tissue volume filled by edema fluid, and open tissue volume holding  PEEP .

Algorithm for parametrization Determination of RD volumes
In a first step, the algorithm determines  0,min and  0,max of each terminal unit depending on the HU type and its location on the pressure-volume curve at PEEP (Figure 5).Normally aerated terminal units are assumed to be in a fully open state, i.e., their expansion from PEEP happens only due to distension and without any recruitment.Therefore, their current reference volume assumed from the CT scan at PEEP,  0,PEEP , equals  0,max (green curve in Figure 5) and can be determined by solving Eq. ( 2) for  tp,PEEP and the known gas volume  PEEP . and  in Eq. ( 2) are input parameters of the algorithm, i.e., presumed quantities that are fit in the outer optimization as described in Section 3.4.2.
The poorly and non-aerated terminal units are (partially) collapsed and assumed to be located somewhere between a fully open and a fully closed state (blue or orange curve in Figure 5).Given this assumption, we can make no direct conclusions about  0,max based on the gas volume in the CT scan (only about  0,PEEP of the terminal units by Eq. ( 2)).We thus use their  tissue as an indication for  0,max that can be reached by a terminal unit in fully open state.Taking the ratios  0,max ∕ tissue of all normally aerated terminal units, we determine the mean and standard deviation of their normal distribution and chose  0,max for the (partially) collapsed terminal units randomly and according to the probabilistic distribution.To complete the RD volume parameters, we further assume for all terminal units  0,min =  coll ⋅  0,max with  coll being a further input parameter to the parametrization algorithm.

Determination of RD time constants
Injured lungs exhibit a change in aeration during expiration that is faster in the dorsal than in the ventral lung when in a supine position 36 .Since this different temporal behavior is very subtle and might also be attributed to the difference in viscoelastic tissue straining due to gravitation, we do not consider pathology or height in the choice of the RD time constants.The time constants  of all terminal units are randomly chosen from a distribution.The type and the parameters of the distribution are input parameters of the algorithm and specifically determined for a patient.

Determination of critical pressures
In addition to the volumes, each terminal unit also requires the specification of the critical pressures  cl,crit,min ,  cl,crit,max , op,crit,min , and  op,crit,max .Given the constant relationship between the critical pressures defined by Δ max−min and Δ op−cl , they are all fixed as soon as the value of one of them is known (Figure 2).Δ max−min and Δ op−cl are both input parameters to the parametrization algorithm.
Before describing the procedure used to determine the critical pressures, we will offer a few basic thoughts: Poorly ventilated terminal units are considered to be partially collapsed at the time of CT recording at PEEP, so they are assumed to take on an unstable state subjected to repetitive intra-tidal RD during normal ventilation.Therefore, the  0,PEEP of a tissue element determined for  PEEP and  tp,PEEP by solving Eq. ( 2), lies somewhere between  0,max and  0,min and offers a clue about the critical pressures.However, due to the time dependence of RD,  0,PEEP is not the actual target volume  0,targ,ee pursued at  tp,PEEP during normal ventilation, but the end-expiratory point of  0 ( , ) reached after the time of expiration  exp during moving towards  0,targ,ee (Figure 6).In a similar manner, during inspiration for the time  insp ,  0,targ,ei is approximated but not fully reached.The  0 ( , ) of a poorly ventilated terminal unit thus continuously moves in transient states during subsequent normal breath cycles and the state at end-expiration is gathered in the CT scan.The initial reference volumes  0,init of a terminal unit at the beginning of inspiration or expiration remain unknown, and it is unclear whether they somehow achieve a steady state.
On the basis of the above consideration and open question, we derive the following mathematical relations in order to eventually determine the critical pressures of poorly ventilated terminal units.
Beginning at end-expiration, we take an initially unknown reference volume  0,init of a terminal unit denoted as  0,1 and assume a constant transpulmonary pressure  tp,endinsp =  endinsp −  pl , which includes the height-dependent gravitational load of a terminal unit, for the time of inspiration,  insp .According to Eq. ( 5), the new reference volume  0,2 at end-inspiration (and  0,init for the following expiration) yields Assumed pressure-and time-dependent steady-state oscillation of the stress-free reference volume of a poorly ventilated terminal unit along normal ventilation cycles.
Again, a constant  tp,PEEP lasting for the time of expiration,  exp , results in a reference volume reading Replacing the exponential expressions with the constants  I =  − insp ∕ insp and  E =  − exp ∕ exp , rearranging the equations, and continuing the stepwise analysis of consecutive inspiration and expiration phases leads to the general expression for even values of  ≥ 4, with  0, denoting the reference volumes at end-inspiration and in the following referred to as  0,insp , and for uneven values of  ≥ 3, with  0, denoting the reference volumes at end-expiration and hereafter referred to as  0,exp .Applying the general relation for geometric series ∑ ∞ =0   = 1 1− , which is always valid for || < 1 and, thus, true in our case ( insp ,  exp > 0 and 0 <  insp ,  exp < ∞ resulting in 0 <  I ,  E < 1), we get Equations ( 11) and ( 12) then converge for  → ∞ to and respectively.Note that  0,1 dropped out in these equations.Therefore, we reach a steady state oscillation between  0,insp and  0,exp that is independent from the initial reference volume assumed at the beginning.Obviously, the condition  → ∞ is not met in reality.However, some iterative calculations demonstrated a good approximation of  0,insp and  0,exp after only a few breathing cycles for realistic mean values of  47,48 .A normal ventilation during CT recording lasting at least several minutes without changes in the ventilator setting should be a period of time long enough to approach this steady state acceptably.
0,exp in Eq. ( 15) should now equal  0,PEEP determined for all poorly aerated terminal units.For this purpose, we use a sampled range for reasonable critical pressures  crit,max and  crit,min and calculate  0,targ,ei and  0,targ,ee for each parameter set according to Eq. ( 3), assuming constant  tp,PEEP and  tp,endinsp during expiration and inspiration, respectively.The parameter set where  0,exp is the closest to  0,PEEP , derived from  PEEP of the terminal unit, gathered in the CT scan, by solving Eq. ( 2), then holds as final  crit,max and  crit,min .
Evaluating also  0,insp helps to define the current path (opening or closing) an element moves on, which eventually sets  crit,max and  crit,min as  op,crit,max and  op,crit,min , or  cl,crit,max and  cl,crit,min , respectively.If  0 reaches  0,insp ≥ ( 1 −  V 0 )  0,max at end-inspiration but not  0,exp ≤ ( 1 +  V 0 )  0,min at end-expiration, then the terminal unit is assigned to move on the closing path, and in all other cases on the opening path.
Regarding the fully open normally ventilated terminal units, we have no actual indications of the relevant critical pressures at hand, especially below which pressure they start to collapse.Thus, apart from ensuring  cl,crit,max <  tp,PEEP , their critical pressures are randomly chosen from a normal distribution with mean  cl and standard deviation  cl , two further input parameters of the global optimization loop.All normally ventilated terminal units are initially moving on the closing path.
The non-aerated terminal units are assumed to have a similar pathological condition to each other and, therefore, similar critical opening pressures.Thus, their critical opening pressures are randomly chosen from a normal distribution with the mean  op and the standard deviation  op -both again input parameters of the outer optimization loop.According to their aeration state at PEEP, we ensure that collapsed terminal units approximate a closed state at end-expiration by setting  op,crit,min >  tp,PEEP .Further, they are not supposed to open entirely during normal ventilation, guaranteed by  op,crit,max >  tp,endinsp .

Optimization of input parameters
As stated, the algorithm described above is based on a set of input parameters which are collected in Table B2.To adjust these quantities, the summed  - tot behavior of all terminal units is matched to the clinical measurements, i.e., to the  - tot points at end-inspiration and end-expiration of a normal breath cycle and to the quasi-static inflation maneuver.For the procedure, we again consider only the elastic component of the terminal units (Eq.( 2)) incorporating pressure-and time-dependent RD.The calibration process thereby also includes the intra-tidal time-and pressure-dependent RD of the terminal units which is continuously present in pathological tissue regions during ventilation 61 , and should therefore be considered.
Specifically, we assume an ideal normal breath cycle at constant  tp,PEEP and  tp,endinsp , specified by the ventilation measurements of the patient and acting on each terminal unit for the duration of  exp and  insp , respectively.We calculate the gas volume  of each terminal unit at end-expiration and end-inspiration which evolves due to the supposed pressure regime and the time, including RD effects.Finally, summing  of all terminal units enables comparison of  tot to the clinical measurement at the two quasi-static points.
This procedure is repeated for the quasi-static inflation maneuver, taking constant transpulmonary pressure courses finely discretized over time between the measured pressure level at the onset ( tp,qs,start ), the peak inspiration ( tp,qs,max ), and the endexpiration ( tp,qs,end =  tp,PEEP ) of the quasi-static inflation maneuver.The time of inspiration  insp,qs and time of expiration  exp,qs of that maneuver is also taken from the clinical measurements.By calculating and summing the  of all terminal units for every small time interval and the corresponding transpulmonary pressure, we produce an ideal pressure-volume curve including the time-and pressure dependent RD of the terminal units, which can be matched to the measured pressure-volume curve of the patient.

Clinical data
As part of generating an exemplary patient-specific model, we used chest CT images and measurements of a critically ill, endotracheally intubated male patient suffering from moderate ARDS.The patient was treated at the operative intensive care unit of the Department of Anesthesiology and Intensive Care Medicine at University Medical Center Schleswig-Holstein, Campus Kiel.All data were provided in an anonymized format.Ethical approval was obtained from the ethics committee of the Medical Faculty in Kiel, and the underlying study was carried out in accordance with the Declaration of Helsinki.Written informed consent was obtained from the patient's legal representative.
The patient chest CT used was recorded at a PEEP of 8 mbar (PEEP8) with 512x512x339 pixels of 0.726562x0.726562x1mm within 24 h before the start of the underlying study.An exemplary axial view of the lung excised from the CT scan is shown in Figure 4.
As part of the study protocol, various ventilation maneuvers were performed including the profiles required for the model calibration described in 3.4.Remaining ventilation maneuvers that did not enter model calibration were used to validate the performance of the model in Section 4.2.

Model generation
We used Mimics and 3-Matic (Materialise, Leuve, Belgium), to segment the patient-specific model geometry and extract dimensions of the trachea, the lobar bronchi, and the lung lobes from the medical images.The space-filling tree growing algorithm yielded the conducting airway tree shown in Figure 3 (left), consisting of 58321 airway segments and having 29161 terminal units attached to the terminal airway segments.

Model calibration
We determined the parameters in the two components of  pl for the patient in MATLAB by non-linear regression.Fitting Eq. ( 6) to the patient's pressure-volume data resulted in the curve of  vol pl depicted in Figure 7.The fitted course of  weight pl (Eq.( 8)) is qualitatively shown in Figure 4.The chosen values for the parameters in  vol pl and  weight pl are listed in Table B1.The input parameters required for the algorithm described in Section 3.4 were optimized manually to match the quasi-static pressure-volume points of the patient.The time constants  of all terminal units were randomly chosen from a quasi-hyperbolic distribution, i.e.,  ∈  unif [0,1] , with the input parameter  , and unif[0, 1] describing uniformly distributed stochastic values between 0 and 1 34 .We adapted the values for Δ max−min and Δ op−cl following the magnitude of variables with similar meaning in other models 38 .Further, the clinically applied and measured pressures entering the model calibration are provided in Table B3.
The  - tot relations resulting from the input parameter values ultimately chosen and provided in Table B2 are shown in Figure 8, together with the related clinical data curves.Note, that for the normal breath cycle the matching is limited to the quasi-static points along the cycle, i.e., to end-inspiration and end-expiration (orange crosses in Figure 8).As a consequence of the assumptions of the parametrization algorithm, there are no resistive effects included in the simplified model underlying the calibration algorithm.Thus, periods of high flow such as in the beginning of inspiration and expiration in a normal breath cycle can not be reproduced adequately.Figure 9 depicts the  op,crit,max of all terminal units (also structured by pathological types), emerging from the algorithm and the given input parameters. op,crit,min ,  cl,crit,min , and  cl,crit,max behave accordingly.

FIGURE 9
The kernel density functions of the  op,crit,max of all terminal units and the subgroups of collapsed (red dashed), poorly-aerated (orange dotted), and normally ventilated (green dashed/dotted) terminal units;  op,crit,min ,  cl,crit,min , and  cl,crit,max behave accordingly.

Simulation protocol
Using the generated model tailored to a patient's lung we simulated 30 minutes of mechanical ventilation of the patient.We applied the airway pressure obtained from the clinical measurements at the airway opening of the simulation model for this purpose.In order to validate our model, we are presenting several time ranges, including various maneuvers, from this simulation study.Note that we are in this case excluding the quasi-static inflation maneuver used for calibration of the model in Section 3 to not simply reproduce the data the model was fed with.Apart from cycles of normal ventilation, the ranges extracted from the study contain the following maneuvers: 1.A quasi-static inflation maneuver that covers a wide pressure range of the patient's respiratory system (Figure 10).At a starting level at zero end-expiratory pressure (ZEEP) and reaching a peak airway pressure of 33 mbar, the maneuver captures a broader pressure range than the quasi-static inflation maneuver used for model calibration, 2. An inspiratory hold maneuver where the phase of inspiration was prolonged by keeping the driving pressure constant for 6 s (Figure 11), allowing an additional volume increase in the lung, inter alia, due to RD; this maneuver is followed by an expiratory hold maneuver, where the flow was kept close to zero, 3. An inspiratory hold maneuver with occluded tracheal airflow, followed by normal ventilation cycles with half the original driving pressure and another inspiratory hold maneuver with flow occlusion at reduced driving pressure (Figure 12), and 4. A decremental PEEP trial, where the PEEP level is gradually reduced from an elevated level, each time started by an inspiratory and ended by an expiratory hold maneuver with flow occlusion (Figure 13).
The time in Figures 10 -13 indicates the chronological sequence of the extracted ventilation periods.halved driving pressure (Figure 12).The same was true of the prolonged periods of expiration, and especially when the airway pressure is decreased toward ZEEP, prior to the quasi-static inflation maneuver (Figure 10).Therefore, the model also captures the volume course in the lower pressure range very well.Note that the clinical volume curve calculated from the measured airflow exhibited an unpredictible baseline drift -potentially due to sensor disturbances caused by the humidified respiratory air -resulting in a disinctively varying end-expiratory volume level.As a result, we compared the simulation and clinical data in a tidal manner only.However, the drift in the quasi-static inflation maneuver was notably strong (Figures 10) when there was a gain in the measured volume of about 0.4 l compared to the volume at PEEP before the maneuver, without causing a similar tendency in the measurements of the actually directly coupled pleural pressure.The close agreement between the clinical and the simulated pleural pressures along the maneuver justifies the validity of our model.Further, at the same pleural pressures during the quasi-static inflation maneuver ( = 203.5 s and  = 250 s) the volume values should also be similar, which is the case for the simulation.
Evaluating the pleural pressure in more detail, the simulation values closely resemble the clinical data.Note that single swallowing actions by the patient caused temporary disturbances in the measurements or failures of the measurement sensor (e.g., see Figures 10 at  = 225 s or 11 at  = 305 s, and 11 at  = 295 s, respectively) 62 .The pleural pressure base level sometimes changed after such events (e.g., starting in Figure 11, enduring in Figure 12, and again reset in Figure 13).The variation in pleural pressure, however, can still serve for comparison despite a varying base level of the pressure 63 , and it does match well for our simulation results and the measurements.RD dynamics The accumulated behavior of alveolar RD can be seen in Figures 10 -13   2), one percent of open tissue corresponds to an additional gas volume of approximately 37 ml at  tp,PEEP (global).We see continuous intratidal RD, and the single maneuvers have remarkable effects on the degree of recruited volume.As expected, the degree of RD in the model is strongly influenced by the level of airway pressure (Figure 13).The effect of time dependence is especially obvious during the longer periods of expiration and inspiration, e.g., prior to the quasi-static inflation maneuver, or during the inspiratory and expiratory hold maneuvers (Figure 11), when there is ongoing recruitment or derecruitment, respectively.
Moreover, stable recruitment of tissue after the quasi-static inflation maneuver is evident by the ∼ 1% increase in the oscillating open tissue compared to the time before the maneuver (Figure 10).However, we see that the model slightly underestimates the tidal volume after the inflation maneuver, i.e., the volume gain is obviously not sufficient.A similar effect was observed in the inspiratory hold maneuvers (Figures 11 and 12) and at the end of the decremental PEEP trial in Figure 13.Hence, although the compliance increases due to the opening of tissue, the overall permanently recruiting volume in the model is insufficient.
Another interesting aspect is that in the decremental PEEP trial the tidal RD exhibited a smaller amplitude in the higher pressure regions than in the lower pressure regions (Figure 13).In addition, reducing the driving pressure apparently also decreases the degree of tidal RD (Figure 12).To give an impression of the regional model behavior for specific ventilation maneuvers, Figures 15 -18 illustrate the local volumetric strains, the current recruitment state (proportion of open reference volume), and the difference in recruitment between specific time points, all for (i) a normal breath (Figure 15), (ii) along the quasi-static inflation maneuver (Figure 16) not used for parametrization, and for periods of (iii) inspiratory and (iv) expiratory hold (Figures 17 and 18, respectively).We calculated the volumetric strain of a terminal unit by   = ( +  tissue )∕( 0,max +  tissue ).Note that the straining of a terminal unit refers to its stress-free state, and not to the functional residual capacity (FRC) or the end-expiratory lung volume (EELV), as it is usually the case in the medical community.We then have the opportunity to evaluate the absolute strain, and not the strain from a specific state.As expected, we see that the heterogeneity in straining is directly linked to the recruitment state of a terminal unit since (partially) closed terminal units generally react in a stiffer manner (see the volumetric strains and proportion of open reference volumes in Figures 15 -18).Note that we obtain a homogeneous strain behaviour across the lung regionally only influenced by the gravitational pressure gradient if we determine the relative strain based on the currently open tissue by   = ( +  tissue )∕( 0 ( , ) +  tissue ), i.e., although the overall strain in a terminal unit may appear low, the open tissue of the terminal unit can experience higher strains.

Local mechanics
As would be expected given the assumptions in model calibration, the tidal recruitment in normal ventilation occurs predominantly in partially collapsed terminal units or at the edges of atelectases, and to a differing degree (Figure 15, bottom).It is further evident that the regional pattern of tidal RD greatly resembles from breath to breath in normal ventilation, in accordance with the mathematical assumptions we made in Section 3.4.
Along the quasi-static inflation maneuver, a high amount of RD occurs spread across the lung (Figure 16), leading to several stably recruiting terminal units, which we previously recognized already in the global RD curves (Figure 10).After significant derecruitment occurring in terminal units all over the lung model when lowering the airway pressure to ZEEP, the model exhibits a high amount of (re-)opening along the inspiratory branch of the maneuver.This even continues very slightly for a few terminal units during expiration from peak pressure back to PEEP (see yellow terminal units in comparison of the proportion of change in reference volume between  = 221.5 s and  = 254.8s in Figure 16).The vast majority of terminal units transitioning during expiration, however, shows a closing behavior (see blue colored terminal units).When evaluating the net impact of the quasistatic inflation between PEEP at the end of the maneuver ( = 254.8s) and the last regular PEEP in normal ventilation before the maneuver ( = 169.9s), the recruitment prevails (see yellow-red colored terminal units in Figure 16, bottom).
Interestingly, we also observe small derecruitment in several terminal units, again close to zero (see cyan colored terminal units in Figure 16, bottom).We traced this phenomenon back to the volume-dependent pleural pressure,  pl : The permanent recruitment of terminal units during the quasi-static inflation increases the total gas volume of the model at PEEP, which again causes a small raise in  pl .The thereby slightly reduced pressure difference across the terminal units leads to some derecruitment, especially in terminal units that are open at PEEP8, but tend to close due to a critical closing pressure close to  tp,PEEP .In total, however, this effect is very small and thus not observable in the global RD curve (Figure 10, bottom).For a more detailed view, the RD between the specified points in time is quantified for all affected terminal units over the lung height in Figure B1 in the appendix.
The effect of time dependence in RD is particularly evident along the inspiratory hold maneuver in Figure 17.At  = 319.7 s, the time when the inspiration of a normal breath cycle would end, the model exhibits a similar RD pattern as in normal ventilation (Figure 17, bottom left, and Figure 15, bottom).The subsequent prolonged inspiration triggers a considerable amount of additional recruitment until  = 319.7 s (Figure 17, bottom right), which appeared already from the global RD curve in Figure 11.
In Figure 18 (bottom left), we see that many terminal units continue to have a higher amount of recruited volume when returning to PEEP at  = 326.9s compared to the end-expiration before the maneuver at  = 318.5 s.Here again, we identify a slight and at first counterintuitive derecruitment in some terminal units similar to the observations along the quasi-static inflation maneuver and again caused by the volume-dependent pleural pressure (for more details see Figure B1).The reverse effect appears at the end of an endured period of expiration at  = 339.1 s (Figure 18, bottom right, and Figure B1).Compared to  = 318.5 s, the majority of terminal units experiences derecruitment at the end of the applied combination of inspiratory and expiratory hold maneuvers, but there are some terminal units slightly opening up.

DISCUSSION AND CONCLUSION
In current clinical practice, ventilation management is performed in a generic manner, inter alia, according to the patient's estimated height and ideal body weight 14 .Such limited attempts to adapt the clinical treatment to a specific patient lack local insight into the patient's lung and disregard the specific pathology, both of which are important aspects that might have a crucial impact on the minimization of VILI, and eventually the reduction of mortality due to ARDS 7,64 .
In this study, we presented a model for lung tissue incorporating a novel approach accounting for the harmful phenomena of OD and repetitive RD, both of which are widely recognized as contributors to the pathogenesis of VILI.The mechanism of RD is implemented by the pressure-and time-dependent variation of stress-free reference volume.By means of the underlying viscoelastic model components, the approach involves time dependence due to both collapse and opening, and tissue resistance.OD is captured by the open reference volume subjected to straining.In principle, overdistension and collapse can occur at the same time in a tissue component.
Concurrently, we proposed a method to apply the newly introduced RD model multi-compartmentally in the framework of an anatomically accurate computational lung model and to tailor the model to specific patients.The algorithmic parametrization of each of the numerous tissue components in the lung model described herein considers the local pathology of the injured organ in the manner indicated by medical imaging.Therefore, a terminal unit reproduces the degree of RD condition of the modelled region as specified by the CT scan, i.e., fully collapsed, poorly aerated (and thus gradually recruited), or normally ventilated at the pressure level of the lung in the recorded image.In order to also capture the pressure-volume behavior of the lung beyond this single, end-expiratory state, the algorithm input parameters are optimized such that the lung model matches the patient's global respiratory mechanics as observed in clinical measurements.The influence of the chest wall on the lung is accounted for by a customized lung volume-dependent pleural pressure boundary condition acting outside the terminal units.
As an application, we deployed the full modeling concept to an example patient suffering from ARDS and receiving mechanical ventilation in the intensive care unit.In simulations of various clinical airway pressure profiles, we tested the model response and compared the resulting tidal volume and pleural pressure to evaluate the predictive capability of the model.Both quantities show very reasonable results and closely follow the measured data (Figures 10 -13).The alveolar RD model elements mimicked repetitive intra-tidal RD, stable recruitment after elevated pressure levels, and transient opening or closing when experiencing changes in the ventilation mode, or along constant pressure profiles.In addition, the parametrization of the exemplary lung model appears plausible when comparing the CT scan to the local model behavior.Moreover, we observe the characteristic of non-normally distributed but scattered critical pressures in the resulting model parameters (Figure 9), which was also observed experimentally 10,35 .
This novel modeling concept demonstrates a variety of promising and significant aspects with regard to future patient-specific treatments: Firstly, the design of the presented alveolar RD model and the method of customized parametrization in a full lung model will enable us to overcome the lack of knowledge and the inaccessibility of biophysical properties of the local lining liquid in air spaces, which emerged in our previous work 37 and has hindered the predictive patient-specific modeling.It is still impossible to obtain the actual physical quantities crucial to pathological behavior in diseased regions of the lung.We are instead using information about the physical properties of lining fluid of deranged lung units intrinsically included in the CT scan by the gray values, e.g., indicating a collapse tendency of a region at the specific pressure state due to the composition of the present fluid.In the proposed model, this information is phenomenologically integrated into the RD parameters.
Secondly -to offer a glimpse into the future potential -, a more realistic reproduction of volumetric behavior and, thus, distribution of air in the lung can enable a more realistic estimation of the surface available for oxygen exchange.However, drawing any conclusions about blood oxygenation will of course require extending the present approach by an adequate model for gas exchange and lung perfusion.
Last but not least, we again stress the ability of the multi-compartment model presented to capture and evaluate two major contributors to VILI, and we further emphasize in this context especially the essential inclusion of time dependence in RD 10 .Individually tailored to a patient's lung geometry, to the overall mechanics of the thoracic cage and to the lung pathologiesbe it the elevated lung weight, or the local collapse tendencies and its dynamics -this model can deliver a measure of the injury potential of ventilation profiles with respect to a specific lung.By changing the reference volume, we can directly access the amount of repetitive reopening which produces harmful shear stresses at airway walls in real life as well as the opening velocity 12,65 , and we can further detect and quantify excessive straining.Previously developed approaches 27,29,31 can serve as estimation for volutrauma and atelectrauma in our model, but in a patient-specific and locally resolved manner.This direction will certainly require further research in multi-patient studies, which is currently ongoing, and a thorough investigation and definition of safe thresholds for the measures mentioned.
To give an even broader perspective on the usability of the novel modeling approach, we outline -just in a nutshell -some potential ideas how such models can enhance mechanical ventilation.One straightforward and evident concept is to leverage such a model to provide additional insight and data concerning non-measurable quantities deep down in the lung during ventilation.Physicians could utilize this information to differentiate between cases and guide subsequent treatment decisions.Additionally, such (so far completely missing) data points could be very helpful when employing novel AI-based approaches in healthcare, as AI relies on the availability of data.The second category of ideas involves utilizing the model to investigate various ventilation maneuvers and settings, enabling the physicians to observe alterations in relevant quantities of interest, and decide for the best option.Moving to the highest level, there exists the possibility to use the model to suggest specific maneuvers or parameters for individual patients or even autonomously optimizing them in the long run.
Finally, we want to address the use of esophageal pressure measurements for the model calibration.While becoming increasingly common, those are not yet a standard practice in the medical management of ventilated patients.In this study, we utilized the measured esophageal pressure due to its inclusion in the example patient's study protocol, and our objective to also reproduce the volume-dependent pleural pressure with the presented lung model.In general, however, it is important to note that model calibration can be realized without direct measurement of the esophageal pressure, especially as long as all pressure-volume relations remain consistent throughout calibration and simulation.For example, in a previous study, we calibrated a simpler model without patient-specific esophageal pressure measurements 26 .However, note that the model then operates on an artificial pleural pressure level, restricting its ability to represent actual lung stresses.As mentioned earlier, measured esophageal pressure qualitatively reflects pleural pressure behavior, although not quantitatively 62,63 .Thus, using measured esophageal pressure might equally introduce a discrepancy in absolute values.Moreover, measurements in the esophagus are strongly influenced by various factors 63 , potentially leading to further inaccuracies.Despite these uncertainties, neither the parametrization nor the simulation results have been significantly compromised in this work.This indicates a certain robustness of the model to variations in pleural pressure.Further, we are often not only or mainly interested in absolute values, but rather in relative quantities, for instance for different ventilation maneuvers.

Limitations
The present model includes a few shortcomings along with the advantages and opportunities described.To begin with, the optimization of input parameters for the described method to calibrate the model is improvable.Regarding the example patient presented, the RD time dependence is obviously not yet properly met by the chosen set of input parameters (e.g., see slope and amount of volume increase along the inspiratory hold maneuver in Figure 11).Therefore, we have to be cautious when interpreting the calibrated values of the critical pressures (Figure 9) and the absolute percentage of tidal RD currently given by the example model.So far, the RD dynamics in the model should only be evaluated qualitatively.We observed that the interaction of single input parameters and their influence on the resulting pressure-volume behavior are to some extent unpredictable.As a solution, the manual adaption of input parameters (Section 3.4.2) can be enhanced, e.g., by an inverse analysis used to automate the model calibration and obtain the best parameter fit by finding a global minimum for the optimization.In this context, we have already presented novel and very promising Bayesian based approaches that may be useful to identify the parameters for the present model 57,58,59 .Doing so enables both pure parameter optimization -maybe even for the full model without using the proposed calibration method to also integrate data from other measuring equipment -as well as considering and testing different distributions of the time constants , e.g., quasi-hyperbolic, exponential, or lognormal, or the pathology-dependent design of variables like the time constants and tissue stiffnesses  26 .The latter goes beyond the scope of this study, but is a valid object of research as remodeling of tissue occurs already in an early stage of ARDS 66 .
Furthermore, while the model introduces RD dynamics, it does not incorporate certain effects caused by the complex fluid mechanics at the air-liquid interface, that are known to appear during RD and that can affect airway walls and cells.These can be seen for example in numerous studies that investigate the complex injury processes leading to atelectrauma and the diverse characteristics that influence it 67 .In particular, the frequency of RD events and the presence of high gradients in wall pressure and shear stress appear to be detrimental 12,65 , with the latter being influenced by various factors.There are interesting approaches in the literature to assess the damage caused by RD 29,68,69 .These and other approaches could provide valuable inspiration for extensions of the proposed model.Our current hypothesis is that these models and underlying insights could be used in an additive manner.One possibility would be to incorporate information on stresses, for instance from moving and rupturing liquid bridges, by adding it to the stress states directly obtained from the proposed model, once the model detects that recruitment occurs.Notably, the presented model might also provide quantities pivotal for determining the extent of RD-induced stress itself, such as opening velocities, or frequency and amount of RD.Naturally, further investigations are necessary to ensure a comprehensive reproduction of the injury potential.
Another limitation in our model is the lack of local mechanical interdependence between the terminal units and of parenchymal tethering between tissue and conducting airway elements.We see an interaction between the terminal units through the coupling by a volume-dependent boundary condition of pleural pressure where the recruitment of some terminal units triggers a slight derecruitment in others, and vice versa.However, the model does not include the direct influence of neighbouring elements onto each other, especially onto their straining behavior when RD occurs.For the airway-tissue interaction, we consider local pressure differences by calculating the external pressure acting on an airway element as the alveolar pressure of the closest terminal unit.We suppose that this already accounts for a large part of the airway-tissue interaction.However, the tethering of the lung parenchyma on the airway wall is neglected.There are promising approaches in literature to model those two types of interdependence 30,50,70 .Though, it is unclear whether the concepts are valid in the current context of RD dynamics and local pathologies, and how to couple them to the model.We therefore refrained from including them in the present work, and point to the need for further investigation of mechanical interdependence as it is assumed to be of importance regarding the phenomenon of RD 25,71 .
Lastly, we purely focused on the alveolar RD model and did not include airway collapse dynamics into the full lung model.Although the approach presented is able to mimic the key features of airway RD, the effect of gas trapping is neglected.That shortcoming can easily be tackled by the additional application of the previously presented airway RD model 37 .However, the problem of a proper parametrization of the critical RD pressures in the airways remains.Therefore, further research on the impact and also the clinical relevance of capturing gas trapping in the model may be advisable.
Despite the aforementioned limitations, the newly introduced alveolar RD components and their integration into a comprehensive multi-compartment model of the lung offer promising and relevant opportunities, as outlined above.In conclusion, the approach presented herein has the potential to individually estimate the benefit or risk of VILI caused by arbitrary ventilation strategies, and to significantly contribute to improving clinical ventilation therapy.

B.2 Calibration of terminal units
The parameters entering the calibration of the terminal units in the lung model of our example patient (Section 3.4) are listed in the following tables.

FIGURE 1
FIGURE 1 Schematic representation of the four-element Maxwell model underlying the terminal units.The non-linear spring  E1 has a variable length which represents the concept of RD by variation of the stress-free reference volume achieved in the present research.See the text herein for a specific definition of all variables and parameters.
,min P cl,crit,min ΔP max-min C lo s in g T U O p e n in g

FIGURE 4
FIGURE 4 Representation of the height-dependent pressure state of terminal units at PEEP, being  tp,PEEP = PEEP −  pl ; the external pressure  pl consists of a volume-dependent component  vol pl and a height-dependent superimposed pressure  g pl

FIGURE 5
FIGURE5 Determination of the location of a terminal unit on the pressure-volume curve and its RD parameters depending on the pathological condition of the terminal unit (indicated by HU) and its individual pressure and volume state at PEEP, i.e., the height-dependent  tp,PEEP and the corresponding gas volume  PEEP extracted from the CT scan.

FIGURE 7
FIGURE 7 Relationship between the volume fraction and the pleural pressure used as boundary condition in the model (light blue) determined from the measured pleural pressure-volume curve (dark blue).

FIGURE 8
FIGURE 8 Idealized pressure-volume behavior (light blue solid) and quasi-static points (orange crosses) of all terminal units of the lung model resulting for the finally chosen input parameters of RD parametrization, compared to clinical measurements of the quasi-static inflation maneuver (dark blue solid) and a normal breath cycle (red dashed) .

FIGURE 10
FIGURE10 Simulation results of the patient-specific computational lung model for a clinically applied airway pressure (top) consisting of normal ventilation and a quasi-static inflation maneuver starting at ZEEP: From top to bottom, the tidal volume and pleural pressure produced by the model (dotted light blue) and measured int the clinic (solid dark blue), and the percentage of open tissue in the whole model (bottom).

FIGURE 11 FIGURE 12
FIGURE11 Simulation results of the patient-specific computational lung model for a clinically applied airway pressure (top) consisting of normal ventilation and an inspiratory hold maneuver without occlusion of the airflow: From top to bottom, the tidal volume and pleural pressure produced by the model (dotted light blue) and measured int the clinic (solid dark blue), and the percentage of open tissue in the whole model (bottom).

FIGURE 13
FIGURE13 Simulation results of the patient-specific computational lung model for a clinically applied airway pressure (top) consisting of normal ventilation with a stepwise decremented PEEP, where each step was initiated by an inspiratory hold maneuver with occluded airflow: From top to bottom, the tidal volume and pleural pressure produced by the model (dotted light blue) and measured int the clinic (solid dark blue), and the percentage of open tissue in the whole model (bottom).
by the percentage of open tissue volume.With the overall open stress-free reference volume of the model of 3.25 l which theoretically gives a gas volume of about 3.69 l at  tp,PEEP (global) when solving Eq. (

Figure 14 FIGURE 14
FIGURE 14 Comparison of the CT scan and the opening proportion of terminal units gradually from fully closed (grey) to fully open (black), in the patient-specific lung model at PEEP ( = 31.6s).
The proportion of open reference volume,  open (), in a terminal unit was determined according to  open () = ( 0 ( , ) −  0,min )∕( 0,max −  0,min ).To evaluate the amount of RD between two states   and  +1 , we computed the differences in the proportion of open reference volume by  open ( +1 ) −  open (  ), such that negative values indicate the closing proportion of  0,max of a terminal unit, and positive values the opening proportion, respectively.Black colored terminal units do not exhibit any change in reference volume between the specified time points.

FIGURE 15
FIGURE 15 Local model behavior during a normal breath cycle at end-expiration ( = 169.9s, 1 ⃝) and end-inspiration ( = 170.8s, 2 ⃝); from top to bottom: volumetric tissue strain in the terminal units, their recruitment state indicated by the proportion of open reference volume  0 (from fully closed = 0 to fully open = 1), and the difference in recruitment between the specified time points (i.e., the proportional change of open reference volume  0 ); view: axial cut through the supine lung model.

FIGURE 16
FIGURE 16 Local model behavior along the quasi-static inflation maneuver at the start ( = 188.2s, 1 ⃝), at the peak ( = 221.5 s, 2 ⃝) and at the end ( = 254.8s, 3 ⃝) of the maneuver; from top to bottom: volumetric tissue strain in the terminal units, their recruitment state indicated by the proportion of open reference volume  0 (from fully closed = 0 to fully open = 1), and the difference in recruitment between the specified time points (i.e., the proportional change of open reference volume  0 ); the difference in recruitment at PEEP8 before and after the maneuver is evaluated between  = 169.6 s and  = 254.8s (bottom); view: axial cut through the supine lung model.

FIGURE 17
FIGURE 17 Local model behavior along the inspiration of an inspiratory hold maneuver at end-expiration before the maneuver onset ( = 318.5 s, 1 ⃝), after the inspiration time of a normal breath cycle ( = 319.7 s, 2 ⃝) and at end-inspiration ( = 324.7 s, 3 ⃝); from top to bottom: volumetric tissue strain in the terminal units, their recruitment state indicated by the proportion of open reference volume  0 (from fully closed = 0 to fully open = 1), and the difference in recruitment between the specified time points (i.e., the proportional change of open reference volume  0 ); view: axial cut through the supine lung model.

FIGURE 18
FIGURE 18 Local model behavior at points of end-expiration before ( = 318.5 s, 1 ⃝) and after ( = 326.9s, 2 ⃝) an inspiratory hold, and at the end of an appended expiratory hold ( = 324.7 s, 3 ⃝); from top to bottom: volumetric tissue strain in the terminal units, their recruitment state indicated by the proportion of open reference volume  0 (from fully closed = 0 to fully open = 1), and the difference in recruitment between the specified time points (i.e., the proportional change of open reference volume  0 ); view: axial cut through the supine lung model.
init and can be visualized by the intersection of the slope of the reference volume at  = 0 (blue dashed line) and  0,targ ( ).See the text and nomenclature for further explanation and a specific definition of all variables and parameters.
τ FIGURE 2 Pressure and time dependent collapse dynamics of the alveolar RD model presented.Pressure dependence: The targeted reference volume  0,targ ( ) of a terminal unit lies between the minimal and maximal reference volumes  0,min and  0,max , respectively, depending on the pressure difference  experienced by the terminal unit, and the active path, i.e., the opening (blue) or closing path (orange) specified by the minimal critical pressures for opening and closing,  cl,crit,min and  op,crit,min , and the pressure offsets Δ max−min and Δ op−cl .Time dependence: Starting at an initial reference volume  0,init ,  0,targ ( ) is approximated over time  based on the time constant , which finally results in the current reference volume  0 ( , ). determines the time delay of the change in  0 ( , ), Table B1 specifies the parameters used for the volume-dependent pleural pressure component  vol pl and the static contribution  weight pl, both introduced in Section 3.3 and fit from clinical data of the example patient.

TABLE B1
Parameters specifying the pleural pressure  pl (Section 3.3) for the example patient.
Table B2 provides the input parameters and Table B3 lists the ventilator setting and ventilation measurements used for the parametrization procedure described.

TABLE B2
Input parameters for RD parametrization (Section 3.4) fitted for the example patient.

TABLE B3
Ventilation parameters of the example patient used for RD parametrization (Section 3.4).FIGURE B1 Proportional change in the reference volume of all terminal units affected by RD between two specified time points; orange squares indicate recruiting and blue circles derecruiting terminal unit, respectively; the coordinates of the lung height range from approximately 70 mm (ventral) to -110 mm (dorsal).Tolerance factor  Stiffness parameter in  E1 equal for all terminal units  cl Mean of  cl,crit,max for normally ventilated terminal units  cl Standard deviation of  cl,crit,max for normally ventilated terminal units  op Standard deviation of  op,crit,min for collapsed terminal units  v ,  v , and  v Patient-specific parameters of  vol pl  w ,  w ,  w ,  w , and  w Patient-specific parameters of  weight pl ℎ balloon Height of the esophageal balloon  pl External pressure depending on  tot and ℎ of a terminal unit Weight-dependent pleural pressure component of  pl on a terminal unit depending on its ℎ  frac Volume share of the gas volume increase in all terminal units  tot,max Maximal gas volume of all terminal units during quasi-static inflation maneuver  tot,PEEP Gas volume of all terminal units at PEEP  cl,crit,max Maximal critical pressure of the closing path of a terminal unit  cl,crit,min Minimal critical pressure of the closing path of a terminal unit  crit,max Maximal critical pressure of a terminal unit  crit,min Minimal critical pressure of a terminal unit  op,crit,max Maximal critical pressure of the opening path of a terminal unit  op,crit,min Minimal critical pressure of the opening path of a terminal unit  Gas flow into a terminal unit  Current gas volume of a terminal unit  0,exp Current reference gas volume of a terminal unit at end-expiration  0,init Initial reference gas volume of a terminal unit  0,insp Current reference gas volume of a terminal unit at end-inspiration  0,max Maximal reference gas volume of a terminal unit  0,min Minimal reference gas volume of a terminal unit  0,PEEP Reference gas volume of a terminal unit at PEEP  0,targ,ee Target reference gas volume of a terminal unit at end-expiration during normal ventilation  0,targ,ei Target reference gas volume of a terminal unit at end-inspiration during normal ventilation  0,targ ( ) Target reference gas volume of a terminal unit  0, Current reference gas volume of a terminal unit at end-inspiration or end-expiration  0 () Stress-free reference gas volume of  of a terminal unit at time   0 ( , ) Current reference gas volume of a terminal unit Transpulmonary pressure of a terminal unit at end-inspiration  tp,PEEP Transpulmonary pressure of a terminal unit at PEEP (applied during the CT recording)  tp,qs,end Transpulmonary pressure of a terminal unit at end-expiration of the quasi-static inflation maneuver  tp,qs,max Transpulmonary pressure of a terminal unit at peak inspiration of the quasi-static inflation maneuver  tp,qs,start Transpulmonary pressure of a terminal unit at the start of the quasi-static inflation maneuver