In this paper, we develop structured tree outflow boundary conditions for modeling the airflow in patient specific human lungs. The utilized structured tree is used to represent the nonimageable vessels beyond the 3D domain. The coupling of the two different scales (1D and 3D) employs a Dirichlet–Neumann approach. The simulations are performed under a variety of conditions such as light breathing and constant flow ventilation (which is characterized by very rapid acceleration and deceleration). All results show that the peripheral vessels significantly impact the pressure, however, the flow is relatively unaffected, reinforcing the fact that the majority of the lung impedance is due to the lower generations rather than the peripheral vessels. Furthermore, simulations of a hypothetical diseased lung (restricted flow in the superior left lobe) under mechanical ventilation show that the mean pressure at the outlets of the 3D domain is about 28% higher. This hypothetical model illustrates potential causes of *volutrauma* in the human lung and furthermore demonstrates how different clinical scenarios can be studied without the need to assume the unknown flow distribution into the downstream region.

## Introduction

Mechanical ventilation is an indispensable tool for the survival of critical care patients suffering from acute respiratory distress syndrome (ARDS) or acute lung injury (ALI). However, the mortality rates of these diseases remains high, approximately 40% (1). In addition, it is widely known that the use of mechanical ventilation is itself the cause of a number of further associated complications (2), which are collectively termed ventilator induced lung injury (VILI). This is despite the more recent adoption of modern *protective* ventilation strategies, where lower tidal volumes are administered while retaining a positive pressure at the end of expiration, which increases lung functional residual capacity and alveolar recruitment. These protocols improve the situation considerably over conventional ventilation (3,4,5). However, there are competing ventilator strategies each with their own merits and drawbacks, for example, the level of positive end-expiratory pressure (PEEP) required by the patient (6,7). Biologically, these aforementioned ventilator related diseases manifest themselves at the alveolar level and are characterized by inflammation of the lung parenchyma, albeit still much controversy surrounds the precise mechanisms (8). VILI can be characterized into separate but not mutually exclusive protocol dependent mechanisms: these are *volutrauma* (alveolar overdistension), *atelectotrauma* (recruitment-derecruitment), and *biotrauma* (inflammation) (9). Predisposition to these conditions is brought about by heterogeneous lung mechanics when the patient is suffering from ARDS or ALI (9), for example, the functional capacity of a diseased lung can be reduced resulting in alveolar overdistension (10). Mortality is usually caused via indirect mechanisms resulting from VILI (11) in which a cascade of events leads to sepsis or multiple organ failure. Understanding the reason why the lungs still become damaged or inflamed during mechanical ventilation is a key question sought by the medical community.

Computational modeling of the tracheobronchial region represents an approach that can provide more in depth knowledge into how the human lung functions mechanically, both during normal breathing and with the assistance of mechanical ventilation. Of particular interest are what conditions within the lung mechanically differ under both diseased and healthy conditions. Currently, a number of limitations exist for modeling the lower airways (in this paper lower airways and higher generations are used synonymously to mean peripheral vessels). First the resolution of CT imaging is restricted (maximum 0.5 mm resolution), hence, smaller vessels are not visible for segmentation. Second, even if the entire lung was fully segmentable, currently, it is not computationally feasible to simulate the full lung tree. Due to these points, alternative methods for modeling the lung must be sought.

Previously, there have been extensive models of nonpatient specific human lungs (12,13,14) based around literature data such as Weibel (15) and Horsfield et al. (16). Furthermore, more recent studies have utilized patient specific CT data (17,18,19) and CT ovine lungs (20) under a variety of boundary conditions, including zero pressure conditions, prescribed flowrate, and impedance. In addition, the simulations of Refs. 17,21 included for the first time the effects of fluid-structure interaction (FSI).

There has been considerable work in the area of lung impedance, specifically related to 1D transmission line models (22,23,24,25). Impedance of the lung is an extremely important phenomenon as it plays a central role in the development and distribution of lung disease. More specifically, the heterogeneity of lung disease is brought about by changes in impedance conditions within the lung. These changes further affect how flow and pressure distribute through the lung. The aforementioned impedance models assume 1D fluid mechanics throughout the entire lung. However, from previous 3D simulations it is well known that the flow in the upper portion of the lower airways is not accurately modeled with 1D equations since typical bifurcating flow phenomena, such as flow recirculation, exist. For this reason, a fully coupled 1D-3D approach must be used and is hence developed in this paper. The developed impedance is based on an asymmetric structured tree methodology as introduced and used for blood in Refs. 26,27,28. This cited work has been modified in order to represent the physiological environment in the pulmonary tree. The coupling of the 1D model to the 3D model is achieved utilizing a Dirichlet to Neumann approach introduced in Ref. 29 for arterial blood flow. The aim of this paper is to apply this modeling approach to lung simulations and also to study the mechanics of the lung under artificial ventilation. In particular, using such an approach, the entire conducting airway can be modeled efficiently obtaining important pressure and flow information.

## Methodology

### 3D Geometry Reconstruction

The geometry used in the simulations was segmented using the commercially available segmentation software MIMICS (Materialise, Leuven, Belgium). The standard CT scans had a resolution of 0.6 mm. This allowed up to seven generations of the bronchial tree to be obtained. Following segmentation the outlets of the geometry were cut normal to the flow direction and the geometry was exported in stereolithography (STL) format. Figure1 shows the lung model with outlets cut normal to the flow direction.

### Meshing

The geometry was then meshed in the commercially available meshing software HARPOON (Sharc, Manchester, UK). However, before meshing, the outlets here extruded within HARPOON in order to satisfy the requirements of the boundary conditions (see Sec. 2). In this software, a high quality mesh in the complex bronchial tree is a matter of minutes, an example mesh with extruded outlet extensions is shown in Fig.2. The final mesh contained approximately $2\xd7106$ tetrahedral elements equating to $1.6\xd7106\u2002degrees$ of freedom.

### Flow Solver

*in-house*multiphysics research code BACI. This code takes advantage of message passing interface (MPI) parallelization, hence, the solution of large scale simulations is achieved in a realistic time frame. The fluid domain was spatially discretised with stabilized finite elements. In time, the equations were discretised using a one-step-$\theta $ scheme. The resulting linear system was solved using the generalized minimal residual iterative solver with multilevel preconditioning, which is implemented in our code using the open-source package Aztec (30). For all simulations a time step size of 0.004s was used and the Newton iterations were considered converged when the residual dropped below $10\u22126$.

### Boundary Conditions

The inflow boundary conditions used in the present simulations are based around normal breathing and mechanical ventilator prescribed flowrate under inspiratory conditions. For normal breathing, a sinusoidal inspiratory profile was utilized, $q(t)=A\u2009sin(\omega t)$, where $A$ represents the flow amplitude and $\omega $ the circular frequency. This was done under two conditions with a mean flow of 0.2 l/s and 0.5 l/s, respectively, (Fig.3a). This basic breathing model in the present study was only used as a benchmark since under nonmechanical ventilation conditions the upper airways (nasal pharynx) should also be included, in particular, because the glottis can produce turbulence in the trachea and the first few generations of the bronchial airways. Under mechanical ventilation conditions, the endotracheal tube is inserted into the level of the trachea, hence, we have not included the upper airways in our simulations.

For mechanical ventilation, a typical protective ventilator constant flow profile, see Fig. 3b, was utilized. This profile is characterized by almost constant flow, apart from a rapid acceleration and deceleration period at the beginning and end of inspiration. Due to the acceleration periods and high frequency oscillations, the mechanical ventilation is expected to lead to different flow and pressure dynamics due to the considerably higher *slew**rates*. The two aforementioned profiles (sinusoidal and constant flow ventilation) cover the situations, which are available to the anaesthetist as control variables for ventilation.

At the walls of the model all velocity components are zeroed in order to enforce the *no slip* condition. Furthermore, since this study wants to focus on the effects of boundary conditions, the walls are assumed to be rigid. Future work will address the issue of FSI as has previously been investigated by our group (17), where simpler boundary conditions were used.

#### Impedance Tree

The above method was implemented into our in-house finite element computational fluid dynamics (CFD) code as part of a outflow boundary condition, which is explained in the following section.

#### 1D-3D Coupling

## Results

### Normal Breathing

The results of the simulations show well documented characteristics for flow in a bifurcating tree, in particular, centripetal acceleration due to curving vessels leads to flow phenomena such as flow recirculation and secondary flow. This all arises from the need to satisfy the continuity condition, leading to Dean like vortices forming in the daughter vessels. Furthermore, the nonuniform geometry causes acceleration of the fluid in the higher generations (for clarification the lowest generation is the trachea), which is not observed in more idealized geometries (see flow vectors in Fig.5).

The narrowing of the vessels in these locations leads to *jetting* of the flow and, hence, pressure decreases. The flow remains relatively disturbed throughout the model, reinforcing the need for 3D simulations in the lower generations of the bronchial tree. In the higher generations, the flow is relatively uniform meaning the application of the impedance boundary conditions in these locations is acceptable.

The most notable change between standard traction free and impedance outflow conditions is the pressure distribution, both spatially and temporally. In Fig.6, the spatial pressure distribution is compared for the two different boundary conditions: traction free (up to seven generation model) and impedance (seven generation model plus up to 13 generations of the 1D tree depending on the outlet size). Evidently, there are marked differences. For the present model, the pressure variation between different outlets did not vary considerably due to the similarity in outlet size. The time history of pressure is completely different for the two different scenarios. From a physiological prospective the above is of utmost importance because determination and understanding of correct airway pressure drop plays an important role in lung mechanics. The above immediately reinforces the need to consider the downstream domain when modeling the tracheobronchial region.

To observe the difference in more detail, normalized values were plotted (Fig.7). This shows that with the impedance conditions there is a slight narrowing of the pressure-flow loop; this indicates an increase in the resistive component of impedance. Overall, the two loops are very similar exhibiting nonlinear behavior with peak pressure and flow basically coinciding. This is to be expected based on the similar outlet sizes throughout the model. The normalized pressure distribution within the 3D domain is effected only marginally by the impedance of the downstream domain. The result is rather interesting, as it suggests that the opposition to flow comes mainly from the upper parts of the airway rather than the impedance of the lower tree. However, the pressure change between impedance and traction free is significant and is likely to play a central role when fluid-structure interaction of the 3D domain is considered.

Figure8 shows the distribution of mean pressure along two different centerline paths: left lobe (gray) and right lobe (black). Overall there is a drop in the pressure descending the vessel. However, due to disturbed flow phenomena, pressure increases are observed in the higher generations. This result is in contrast to more idealized models, which lack the pressure variation due to local curvature changes. Hence the pressure in such models tends to drop in a relatively uniform manner along the vessel. The results here also indicate the strong dependence on asymmetry, which previously has been reported to be an important factor in the resulting pressure distribution. In particular the distribution between the left and right lung lobe bear no relationship. Furthermore, when the results are compared with the traction free model (data not shown) greater spatial gradients in pressure are observed. Finally, it is seen from the results that the largest pressure drop is observed in the higher generations.

The effect of different trees on the pressure distribution in the 3D domain is significant, however, this is to be expected as the impedance is dependent on the size of the tree, i.e., the smaller the tree the greater the impedance. This is demonstrated in Fig.9, where two different symmetric Weibel trees ($\alpha =\beta =0.5$ and 0.9) are attached to the end of the model. Evidently the pressure in the 3D domain is considerable higher. In particular the resulting pressure in the faster diminishing tree is higher because essentially the tree can accommodate less flow (due to reduced cross-sectional area). This represents a constriction of the entire lung tree.

The above suggests that the depth of the tree is the most important factor, i.e., the number of generations. This means the structure of the downstream tree is *relatively* insignificant to the spatial pressure distribution, especially the difference between symmetric and asymmetric trees. Overall, the models utilizing an asymmetric tree are better because they are a closer representation of in vitro measured trees, however, some experimental validation based around pressure measurements is still required. In all presented models the effects on flow appear to be relatively minimal, there is a small change in the field, however, this appears to be relatively insignificant.

Finally, we consider how the boundary conditions can be used to simulate *hypothetical* disease in the human lung. The current impedance model has the advantage that no a priori assumptions need to be made about the flow since the flow in the 3D domain is governed by the peripheral vessels. This is a major advantage over previous impedance implementations. Here, we consider a light constriction of the left superior lobe, i.e., a lung of reduced functional residual capacity. This could represent clinical scenarios such as alectasis or a derecruited region of the lung. This model is only indicative and is not meant to reproduce the complex nature of true lung disease since lung disease formation tends to be heterogeneous. The constriction of the airways is modeled by increasing the impedance on the outlets. Here, the model responds by letting less flow into this lobe and naturally diverting flow to other regions of the lung. This shows that with our method we can simulate heterogeneous changes in the lung. Furthermore, considering Fig.10, the pressure-flow relationship is only really affected in the vicinity of the occlusion. In Figs. 10a,10b, it is evident that there is a widening of the loop (hysteresis) indicating an increase in the phase difference between pressure and flow. However, in other parts of the lung the relationship remains very similar (between disease and undiseased scenarios). This suggests the response is quite local, i.e., pressure and flow is increased throughout the lung, however, the dynamic between the two remains similar.

### Mechanical Ventilation

Under mechanical ventilation complex flow dynamics result from both geometric and temporal influences. Ultimately geometry is the major factor in the resulting flow distribution, however, the rapid deceleration observed late in the inspiratory phase results in larger recirculation zones forming in the higher generations. The flowrate difference between the two different boundary conditions shows very similar profiles (data not shown). In addition, the instantaneous flowrate difference is maximal 4% higher under impedance conditions. The most notable change between standard traction free and impedance outflow conditions is again the pressure distribution, both temporally and spatially. In particular the max pressure difference between the two different boundary conditions at peak inspiration is 27% in the second generation (traction free condition is 27% lower than impedance). The observed temporal gradients in pressure are interestingly high in the higher generations, however, the dynamics tended to be similar. This effect is driven by the nonuniform geometry.

The actual maximum pressure level for the first four to five generations generally varies only a small amount. Interestingly in the left lobe of the lung, the pressure over the cycle in the third generation is higher than in the second. This is driven by the flow route, more air is diverted into the left superior lobe, where the pressure is reduced (data not shown). This results in less air flowing into the left inferior lobe. The flow route is primarily driven by the curvature of the second generation of the left hand side of the lung, where the air is skewed toward the interior surface of the lung, hence flow continuity results in secondary currents, which lead to flow diverting back toward the anterior surface. This aforementioned phenomenon is primarily due to the nonuniformity of the geometry, meaning the cross-sectional area change plays an important role in the pressure and flow distribution. This particularly suggests that models based around circular cross sections fail to capture essential dynamics such as pressure variations due to cross-sectional changes.

Figure11 compares the pressure-flow relationship between ventilation and breathing. Clearly there is a change in this relationship. This is evidenced by the pressure-flow loop area increase, hence indicating that inertial effects are important (pressure is leading flow). This relationship is also seen to be much more nonlinear. This is specifically true during the acceleration phases. The small lag that does occur comes about due to transient fluid inertia, which results directly from higher slew rates observed during this ventilation protocol. Ultimately the flow in the lung can be characterized as *quasi-steady* (Womersley numbers $\u223c1$), which means the flow responds to the changing pressure gradient in a very rapid fashion. Both boundary conditions, under mechanical ventilation, exhibit phase differences; this indicates that the upper generations contribute more dominantly to the lung impedance than the peripheral vessels, this is also backed up by the dynamics being not significantly different. However, the peripheral vessels are obviously important for correct pressure distribution/level in the lower generations, as indicated by the large difference between pressure levels.

Again a *hypothetical* diseased lung is simulated. The results are shown in Fig.12. Most interesting is the result in the superior left lobe (Fig. 12b), where the increased impedance on the branch results in an in phase relationship (between pressure and flow) during acceleration and deceleration. This particularly shows a very different relationship to that of the healthy profile and highlights that the method of ventilation is important.

Modeling this reduced functional lung results in an overall increase in pressure in the lung, see Fig.13. In particular, the mean pressure plotted along a centerline (gray centerline from Fig. 8) is elevated throughout and at the outlet (located in the inferior left lobe) it is approximately 28% higher. The pressure is up to 38% higher in the vicinity of the occluded area. Similar observations are also seen in the right lobe. These results mean that the nondiseased part of the ventilated lung is potentially at risk from *volutrauma* due to an increase in lobular pressure over a normal functioning lung.

## Conclusion

In this paper we have developed and used structured tree outflow boundary conditions for modeling the lower airways of the lung tree up to 17 generations. The coupling of this 1D model to the 3D domain uses a Dirichlet to Neumann approach. This approach is most general as it makes no assumption about the form of the downstream domain. Furthermore, the method is applied in a mathematically well posed formulation, i.e., avoiding a full Dirichlet problem. Currently the impedance model is based purely around a structured tree, however, the present implementation allows essentially any tree to be attached to the end of the 3D domain, for example, a tree developed from a space filling algorithm (37).

First the boundary conditions were compared under normal breathing conditions, demonstrating that there is a clear difference between the pressure levels. In addition, different trees were attached to the outlets, specifically, a Weibel tree with different scaling coefficients. This showed that the pressure level can change significantly. The faster scaling tree coefficient $(\alpha =\beta =0.5)$ represents a tree of reduced volume, causing elevated pressure in the lungs. This highlights the importance of the small airways on the pressure levels in the lower generations. Utilizing models in this way allows the effects of small airways to be investigated, which otherwise remains difficult (38).

We further investigated the effects of constant flow ventilation on lung mechanics. The high slew rate ventilation curve resulted in larger recirculation zones in higher generations, particularly during the deceleration phase. The higher frequencies of the mechanical ventilation profile led to a different pressure-flow dynamic, in particular during the rapid acceleration/deceleration phase and where the flow became slightly negative (compare Fig. 11). The wider pressure-flow loop indicates a phase lag between pressure and flow; this is in contrast to normal breathing, where the pressure and flow remain essentially in phase. This different dynamic is due to the effects of transient fluid inertia. However, *generally* the flow in the lungs is quasi-steady since there is only a small lag between pressure and flow, i.e., flow responds rapidly to the changing pressure gradient. This is to be expected given the Womersley number is only marginally above one for the lower generations and below one for higher generations.

Overall, for all models the flowrate did not change considerably. The fact that the impedance does not change the flow magnitude significantly suggests that models, which are interested in flow related phenomena only, such as particle distribution, potentially do not require such complicated impedance type boundary conditions. At the very least, such simulations can consider only resistance at the outlet terminals. However, in making such an assumption, a significant number of generations must be considered in the model in order to capture the relevant pressure-flow dynamics, i.e., the impedance is predominantly from the upstream area of the model. This is supported by previous medical observations, where it is reported that around 10% of the resistance comes from the small airways (39,40). From a flow physics perspective, this does not mean the number of generation levels, rather it requires that the Reynolds number at the outlet is sufficiently low.

Modeling diseased lungs is essential for understanding ventilation of diseased patients and can also provide insight into how different mechanical maneuvers affect the diseased lungs. Understanding how the lungs respond under mechanical ventilation is currently an unresolved question still trying to be elucidated by the medical community (41). Here, a partial obstruction of the superior left lobe was investigated in order to see the effects on the rest of the lung. This simulation also demonstrated the usefulness of the coupling method, where no assumption is made a priori about the flow into the peripheral airways, i.e., the flow is realistically governed by the impedance in this downstream region. The obstruction was not meant to model the complexity of lung disease but instead see the effects of the peripheral vessels on the lung dynamics. However, it could be indicative of a downstream derecruited region. The results demonstrated for both normal breathing and mechanical ventilation that the major changes in pressure-flow dynamics tended to be localized to the vicinity of the region, where the constricted vessels are located, i.e., the other regions of the lung showed increased pressure and flow, however, the relationship between the variables is virtually unchanged. In these other regions, the mean outlet pressure was on average approximately 28% higher, which could potentially result in overdistension (volutrauma) of the downstream alveoli due to elevated lobular pressure. The model here, although completely idealized, demonstrates that diseased lungs respond differently to healthy lungs. Further work is required in this area to understand heterogeneous lung disease. This model represents a step in the direction to understand the complex lung mechanics occurring during mechanical ventilation. One of the major areas, which requires understanding is the effect of different mechanical maneuvers on lung mechanics in a diseased state.

In literature, alternative approaches related to modeling the total tracheobronchial airways have been proposed (42,43). In Ref. 42, a 15 generation model of the airway system was simulated by calculating the flow in a 3D Weibel triple bifurcation and then repeatedly reproducing this section in parallel, where the outflow from each section becomes the inflow for the next downstream bifurcation unit. In Ref. 43, an idealized 17 generation model based on an anatomical data set was simulated under steady state conditions, however, their mesh used in the model was very course. Although these two methods provide some useful insights into flow structures in the conducting airways our approach has a number of advantages: We can model the whole conducting airway utilizing a realistic geometry for the first seven generations and then attach 1D trees as boundary conditions onto the 3D domain; the 1D domain takes into account the effects of wall movement; we are utilizing transient simulations; and we have a framework in which we have demonstrated that we can model diseased state conditions.

Future work will also encompass the expiration phase, which is nontrivial from a computational perspective. The simulations will also take into account FSI effects using our vast experience in FSI and the efficient implementations in our in-house code (21). This is expected, given the change observed in pressure under impedance conditions to lead to different flow and stress dynamics in the tracheobronchial region. In particular significant changes are expected to be seen in the higher generations due to the reduced wall cartilage content. Currently the 1D model is attached to the outlet of the 3D model obtained from the CT data, future models will include artificial bifurcations based on morphological data to extend the 3D model a few extra generations, this is in order to limit any error associated with using a 1D model. Ultimately, we are working toward a whole lung model, hence, work is currently being undertaken to couple the model presented here to models of the alveolar region (44,34), therefore we will have a model of the tracheobronchial tree coupled to the acinar region. This is expected to alter the results since this is a major part of the respiratory system. This is in the hope that further understanding of the lung, specifically the alveoli under healthy and diseased conditions can be ascertained. Finally, the present models have provided detailed insight into the effects of the peripheral airways on breathing mechanics and also the use of coupling for lung simulations.

## Acknowledgment

We acknowledge the support of Deutsche Forschungsgemeinschaft (DFG) through Project No. WA-1521/8. The CT images of the tracheobronchial tree have been provided by Dr. med Ulrich Barkow and Dr. med. univ. Andreas Strohmaier from “Diagnostische Radiologie,” Stuttgart, Germany.