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.

Figure 1
Segmented lung model (up to seven generations) used in the numerical simulations
Figure 1
Segmented lung model (up to seven generations) used in the numerical simulations
Close modal

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×106 tetrahedral elements equating to 1.6×106degrees of freedom.

Figure 2
Geometry with computational mesh (a) enlarged portion of the trachea and bronchi and (b) higher generations also showing outlet extensions required for better modeling of the outflow boundary conditions
Figure 2
Geometry with computational mesh (a) enlarged portion of the trachea and bronchi and (b) higher generations also showing outlet extensions required for better modeling of the outflow boundary conditions
Close modal

Flow Solver

The governing equations for airflow in the tracheobronchial region are the time dependent incompressible Navier–Stokes equations.
1
2
where u represents the velocity vector (u=[u,v,w]), p is the air pressure, ρ is the density, and μ is the dynamic viscosity. The numerical simulations were performed in our 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-θ 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 106.

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)=Asin(ωt), where A represents the flow amplitude and ω 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.

Figure 3
Inflow curves utilized in the simulations (a) light breathing (Qmean=0.2 l/s) and light activity (Qmean=0.5 l/s) and (b) constant flow mechanical ventilation. Note the ventilation curve was chosen based on the high slew rates and oscillations, which occur in the profile to understand how these affect lung mechanics.
Figure 3
Inflow curves utilized in the simulations (a) light breathing (Qmean=0.2 l/s) and light activity (Qmean=0.5 l/s) and (b) constant flow mechanical ventilation. Note the ventilation curve was chosen based on the high slew rates and oscillations, which occur in the profile to understand how these affect lung mechanics.
Close modal

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 slewrates. 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 airways beyond the 3D domain are modeled using a 1D approach previously described for arterial flow (26,27). This 1D approach is modified in order to model the pulmonary tree rather than the peripheral arterial tree. In order to calculate the impedance of the pulmonary tree, the airways are modeled as 1D flexible elastic tubes. To achieve this, the linearized 1D fluid momentum equation is used and is given by
3
where u is the velocity in the streamwise direction, t time, ρ density, ν viscosity, and r the radius of the vessel. The approximation of linearized 1D momentum is reasonable for the pulmonary tree due to the high wave speed. In addition to the momentum equation, we also have the 1D continuity equation
4
where p is the pressure, q is the flowrate, and C is the vessel compliance, which is given by the following linearized state equation (22).
5
with r as the vessel radius, E as the elastic modulus, and h as the wall thickness. The wall thickness values are determined, based on vessel diameter, from the human Horsfield lung model given in Ref. 24.
Equations 3,4 can then be formulated into an expression for impedance in the frequency domain by utilizing Womersley’s solution, which is an analytical solution for oscillatory flow in a circular tube (31). The derived impedance equation expresses the impedance at the root of each airway segment (Z(0,ω)) as a function of the impedance in the downstream airway tree (Z(L,ω)) and is given by
6
where g is the product of the wave speed c and the vessel compliance C, L is the vessel length, and ω the circular frequency. The wave speed c is given by
7
where A0 is the cross-sectional area, Fj is a function of the zeroth and first order Bessel functions as calculated from Womersley’s solution, and ρ is the density and compliance.
The impedance of each individual branch is then summed in series and parallel in order to obtain the total impedance of a specific tree. At bifurcation points, the impedance of the parent vessel is related to the daughter vessels using a standard bifurcation condition.
8
where subscripts P, L, and R represent parent, left daughter, and right daughter vessels, respectively. In addition, scaling of the tree is achieved by introducing a relationship between parent and daughter vessels.
9
where α and β represent parameters for tree asymmetry. For example, if α=β=constant, the tree would represent a Weibel lung tree (15). In this paper, both symmetric and asymmetric trees were investigated. Asymmetric branching was based on averaged reported values in literature (32,33) with α=0.86 and β=0.7. Recursively applying Eqs. 6,8 from the root of the tree down to a critical radius provides the impedance of the lung tree in the frequency domain. In the current model, the impedance at the root of the tree is included but in all numerical examples in this study it was considered to be zero. Realistic nonzero impedance at the root of the tree will be addressed in future models, hence, taking into account peripheral compliance. This will be achieved by coupling to coupled multiscale simulations of the alveolar region (34 35).
To calculate the time varying pressure at the root of a specific tree, the impedance must be first transformed into the time domain via
10
Finally, a time varying pressure is obtained at the root of the artificial tree by convoluting the impedance z(t) with the history of flow at a specific outlet, i.e.,
11

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

The coupling of the impedance model to the 3D domain was achieved using a previously described Dirichlet–Neumann approach, see Ref. 29 for full details. This method has the advantage of being mathematically well posed and has proven to be numerically stable. In addition, it also has the advantage that the solution of 1D and 3D domain is achieved simultaneously. The method was implemented both explicitly and implicitly into our finite element CFD code, however, little difference was observed between the methods and the explicit version was utilized due to the faster solution times. In brief, given the downstream domain pressure, p(t) (obtained from the impedance model), the time dependent pressure at the outflow of the 3D domain can be represented by
12
where Γs is the outflow boundary surface, w is a weighting function, Ĩ is the identity tensor, τ̃ is the viscous stress tensor, and n is the outward normal. The second term (τ̃) is assumed to be zero, which inherently assumes that the flow is fully developed (τ̃=0). To satisfy this condition, the outlets of the 3D model are extruded. All outlet extrusions were deemed suitably long such that the observed profile was parabolic and streamwise independent. Furthermore, the effect that these extensions have on the pressure is minimal. This was further verified using the traction free condition with and without the extensions and comparing the pressure at the equivalent location. Finally, the use of these pressure based boundary conditions, in particular, how the flow distributes, was verified by looking at the mean flow division ratio as a percentage of the inlet flow. The following distribution was obtained: right lobe (57%), left lobe (42%), right superior lobe (19%), right middle lobe (10%), right inferior lobe (28%), left superior lobe (24%), and left inferior lobe (18%). These flow divisions were deemed suitable based on reported literature values (16). The major advantage of this method is that no assumption needs to be made about the flow, hence, flow distribution is completely governed by the downstream region. Further validation was provided by considering the pressure drop calculated with this coupled approach, we considered the pressure drop over the whole model for the varying flowrates. This instantaneous drop was approximately 33 Pa for a flowrate of 0.5 l/s; if we consider the equivalent values reported in literature (36), the value shows good agreement. The bronchial pressure drop over the acceleratory phase is given in Fig.4.
Figure 4
Pressure flow relationship at the inlet to the trachea during the acceleratory part of the inspiratory phase (normal breathing). The pressure drops are similar to the values reported in literature.
Figure 4
Pressure flow relationship at the inlet to the trachea during the acceleratory part of the inspiratory phase (normal breathing). The pressure drops are similar to the values reported in literature.
Close modal

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).

Figure 5
Zoomed in area of the bronchial airway demonstrating flow vectors at max inspiration for light activity. The flow is characterized by complex recirculation and skewing of high velocity regions to inner walls of bifurcations due to the vessel branching.
Figure 5
Zoomed in area of the bronchial airway demonstrating flow vectors at max inspiration for light activity. The flow is characterized by complex recirculation and skewing of high velocity regions to inner walls of bifurcations due to the vessel branching.
Close modal

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.

Figure 6
Comparison of pressure contours, for light activity, at maximum inspiration: (a) traction free boundary condition (up to seven generations) and (b) impedance boundary condition(up to 17 generations). With the impedance conditions the maximum pressure is approximately 44% higher. Note the lower limit of the scale is set to 5 Pa for visual comparison purposes and for the traction free condition the pressure at the outlets is virtually zero, meaning the pressure in the rest of the tree would be negative.
Figure 6
Comparison of pressure contours, for light activity, at maximum inspiration: (a) traction free boundary condition (up to seven generations) and (b) impedance boundary condition(up to 17 generations). With the impedance conditions the maximum pressure is approximately 44% higher. Note the lower limit of the scale is set to 5 Pa for visual comparison purposes and for the traction free condition the pressure at the outlets is virtually zero, meaning the pressure in the rest of the tree would be negative.
Close modal

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.

Figure 7
Normalized pressure-flow relationship in the bronchial airway. Evidently the pressure and flow are consistently in phase throughout the cycle. This is more noticeable for impedance conditions, where a narrowing of the loop is observed. For the impedance condition, the temporal development of pressure is more rapid. The difference in flowrate between the two conditions is marginal with the flowrate being overall reduced under impedance conditions (data not directly shown, however, at each point along the two above curves the flow is nearly identical).
Figure 7
Normalized pressure-flow relationship in the bronchial airway. Evidently the pressure and flow are consistently in phase throughout the cycle. This is more noticeable for impedance conditions, where a narrowing of the loop is observed. For the impedance condition, the temporal development of pressure is more rapid. The difference in flowrate between the two conditions is marginal with the flowrate being overall reduced under impedance conditions (data not directly shown, however, at each point along the two above curves the flow is nearly identical).
Close modal

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.

Figure 8
Mean pressure for light activity along two different centerlines leading into the left (gray) and right (black) lobes of the lung
Figure 8
Mean pressure for light activity along two different centerlines leading into the left (gray) and right (black) lobes of the lung
Close modal

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 (α=β=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.

Figure 9
Pressure distribution for light breathing plotted along a centerline of the right lobe for two different diameter diminishing ratios (α=β=0.5 and α=β=0.9) in a Weibel tree. The pressure is higher in the faster diminishing tree due to higher impedance levels calculated in the downstream domain.
Figure 9
Pressure distribution for light breathing plotted along a centerline of the right lobe for two different diameter diminishing ratios (α=β=0.5 and α=β=0.9) in a Weibel tree. The pressure is higher in the faster diminishing tree due to higher impedance levels calculated in the downstream domain.
Close modal

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.

Figure 10
Pressure flow distribution for light activity in different generations of the lung under healthy and diseased conditions (partial occlusion of the superior left lobe): (a) second generation left lobe, (b) third generation superior left lobe, (c) third generation inferior left lobe, and (d) third generation inferior right lobe.
Figure 10
Pressure flow distribution for light activity in different generations of the lung under healthy and diseased conditions (partial occlusion of the superior left lobe): (a) second generation left lobe, (b) third generation superior left lobe, (c) third generation inferior left lobe, and (d) third generation inferior right lobe.
Close modal

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 1), 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.

Figure 11
Comparison of normalized pressure versus flow relationship for normal breathing and mechanical ventilation. Evidently the higher frequency ventilation results in a widening of the pressure-flow loop implying a greater phase difference between pressure and flow.
Figure 11
Comparison of normalized pressure versus flow relationship for normal breathing and mechanical ventilation. Evidently the higher frequency ventilation results in a widening of the pressure-flow loop implying a greater phase difference between pressure and flow.
Close modal

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.

Figure 12
Pressure flow distribution in different generations of the lung under healthy and diseased conditions (partial occlusion of the left superior lobe): (a) second generation left lobe, (b) third generation superior left lobe, (c) third generation inferior left lobe, and (d) third generation inferior right lobe
Figure 12
Pressure flow distribution in different generations of the lung under healthy and diseased conditions (partial occlusion of the left superior lobe): (a) second generation left lobe, (b) third generation superior left lobe, (c) third generation inferior left lobe, and (d) third generation inferior right lobe
Close modal

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.

Figure 13
Comparison of pressure distribution (t=0.5 s) under healthy and diseased conditions (superior left lobe occluded). The pressure distribution is elevated throughout the model, hence, overdistension of the alveoli is potentially possible.
Figure 13
Comparison of pressure distribution (t=0.5 s) under healthy and diseased conditions (superior left lobe occluded). The pressure distribution is elevated throughout the model, hence, overdistension of the alveoli is potentially possible.
Close modal

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 (α=β=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.

1.
Rubenfeld
,
G. D.
,
Caldwell
,
E.
,
Peabody
,
E.
,
Weaver
,
J.
,
Martin
,
D. P.
,
Neff
,
M.
,
Stern
,
E. J.
, and
Hudson
,
L. D.
, 2005, “
Incidence and Outcomes of Acute Lung Injury
,”
N. Engl. J. Med.
0028-4793,
353
, pp.
1685
1693
.
2.
Ranieri
,
V. M.
,
Suter
,
P. M.
,
Tortorella
,
C.
,
Tullio
,
R. D.
,
Dayer
,
J. M.
,
Brienza
,
A.
,
Bruno
,
F.
, and
Slutsky
,
A. S.
, 1999, “
Effect of Mechanical Ventilation on Inflammatory Mediators in Patients With Acute Respiratory Distress Syndrome: A Randomized Controlled Trial
,”
JAMA, J. Am. Med. Assoc.
0098-7484,
282
, pp.
54
61
.
3.
Brower
,
R. G.
, and
Fessler
,
H. E.
, 2000, “
Mechanical Ventilation in Acute Lung Injury and Acute Respiratory Distress Syndrome
,”
Clin. Chest Med.
0272-5231,
21
(
3
), pp.
491
510
.
4.
Amato
,
M. B.
,
Barbas
,
C. S.
,
Medeiros
,
D. M.
,
Magaldi
,
R. B.
,
Schettino
,
G. P.
,
Lorenzi-Filho
,
G.
,
Kairalla
,
R. A.
,
Deheinzelin
,
D.
,
Munoz
,
C.
,
Oliveira
,
R.
,
Takagaki
,
T. Y.
, and
Carvalho
,
C. R.
, 1998, “
Effect of a Protective-Ventilation Strategy on Mortality in the Acute Respiratory Distress Syndrome
,”
N. Engl. J. Med.
0028-4793,
338
, pp.
347
354
.
5.
Hickling
,
K. G.
,
Henderson
,
S. J.
, and
Jackson
,
R.
, 1990, “
Low Mortality Associated With Volume Pressure Limited Ventilation With Permissive Hypercapnia in Severe Adult Respiratory Distress Syndrome
,”
Intensive Care Med.
0342-4642,
16
, pp.
372
377
.
6.
Brower
,
R. G.
,
Lanken
,
P. N.
,
MacIntyre
,
N.
,
Matthay
,
M. A.
,
Morris
,
A.
,
Ancukiewicz
,
M.
,
Schoenfeld
,
D.
, and
Thompson
,
B. T.
, 2004, “
Higher Versus Lower Positive End-Expiratory Pressures in Patients With the Acute Respiratory Distress Syndrome
,”
N. Engl. J. Med.
0028-4793,
351
(
4
), pp.
327
336
.
7.
Michaud
,
G.
, and
Cardinal
,
P.
, 2003, “
Mechanisms of Ventilator-Induced Lung Injury: The Clinician’s Perspective
,”
Crit. Care
1364-8535,
7
, pp.
209
210
.
8.
MacCallum
,
N.
, and
Evans
,
T.
, 2004, “
Acute Lung Injury
,”
Anaesth. Intensive Care
0310-057X,
5
(
11
), pp.
389
391
.
9.
Moloney
,
E. D.
, and
Griffiths
,
M. J. D.
, 2004, “
Protective Ventilation of Patients With Acute Respiratory Distress Syndrome
,”
Br. J. Anaesth.
0007-0912,
92
(
2
), pp.
261
270
.
10.
Vieira
,
S. R.
,
Puybasset
,
L.
,
Richecoeur
,
J.
,
Lu
,
Q.
,
Cluzel
,
P.
,
Gusman
,
P. B.
,
Coriat
,
P.
, and
Rouby
,
J.
, 1998, “
A Lung Computed Tomographic Assessment of Positive End-Expiratory Pressure-Induced Lung Overdistension
,”
Am. J. Respir. Crit. Care Med.
1073-449X,
158
, pp.
1571
1577
.
11.
Ware
,
L. B.
, and
Matthay
,
M. A.
, 2000, “
The Acute Respiratory Distress Syndrome
,”
N. Engl. J. Med.
0028-4793,
233
(
2
), pp.
309
319
.
12.
Green
,
A. S.
, 2004, “
Modelling of Peak-Flow Wall Shear Stress in Major Airways of the Lung
,”
J. Biomech.
0021-9290,
37
(
5
), pp.
661
667
.
13.
Liu
,
Y.
,
So
,
R. M. C.
, and
Zhang
,
C. H.
, 2002, “
Modeling the Bifurcating Flow in a Human Lung Airway
,”
J. Biomech.
0021-9290,
35
, pp.
465
473
.
14.
Zhang
,
Z.
, and
Kleinstreuer
,
C.
, 2004, “
Airflow Structures and Nano-Particle Deposition in a Human Upper Airway Model
,”
J. Comput. Phys.
0021-9991,
198
, pp.
178
210
.
15.
Weibel
,
E. R.
, 1963,
Morphometry of the Human Lung
,
Springer
,
New York
.
16.
Horsfield
,
K.
,
Dart
,
G.
,
Olson
,
D. E.
,
Filley
,
G. F.
, and
Cumming
,
G.
, 1971, “
Models of the Human Bronchial Tree
,”
J. Appl. Physiol.
8750-7587,
31
, pp.
207
217
.
17.
Wall
,
W. A.
, and
Rabczuk
,
T.
, 2008, “
Fluid Structure Interaction in Lower Airways of CT-Based Lung Geometries
,”
Int. J. Numer. Methods Fluids
0271-2091,
57
, pp.
653
675
.
18.
Baoshun
,
M.
, and
Lutchen
,
K. R.
, 2006, “
An Anatomically Based Hybrid Computational Model of the Human Lung and Its Application to Low Frequency Oscillatory Mechanics
,”
Ann. Biomed. Eng.
0090-6964,
14
, pp.
1691
1704
.
19.
Lin
,
C. L.
,
Tawhai
,
M. H.
,
McLennan
,
G.
, and
Hoffman
,
E. A.
, 2007, “
Characteristics of the Turbulent Laryngeal Jet and Its Effect on Airflow in the Human Intra-Thoracic Airways
,”
Respir. Physiol. Neurbiol.
1569-9048,
157
, pp.
295
309
.
20.
Kabilan
,
S.
,
Lin
,
C.
, and
Hoffman
,
E. A.
, 2006, “
Characteristics of Airflow in a CT-Based Ovine Lung: A Numerical Study
,”
J. Appl. Physiol.
8750-7587,
102
, pp.
1469
1482
.
21.
Küttler
,
U.
,
Gee
,
M.
,
Förster
,
C.
,
Comerford
,
A.
, and
Wall
,
W. A.
, 2010, “
Coupling Strategies for Biomedical Fluid-Structure Interaction Problems
,”
International Journal for Numerical Methods in Biomedical Engineering
,
26
, pp.
305
321
.
22.
Suki
,
B.
,
Habib
,
R. H.
, and
Jackson
,
A. C.
, 1993, “
Wave Propagation, Input Impedance, and Wall Mechanics of the Calf Trachea From 16 to 1,600 Hz
,”
Am. J. Physiol.
0002-9513,
75
, pp.
2755
2766
.
23.
Lutchen
,
K. R.
,
Hantos
,
Z.
,
Petak
,
F.
,
Adamicza
,
A.
, and
Suki
,
B.
, 1996, “
Airway Inhomogeneities Contribute to Apparent Lung Tissue Mechanics During Constriction
,”
Am. J. Physiol.
0002-9513,
80
, pp.
1841
1849
.
24.
Gillis
,
H. L.
, and
Lutchen
,
K. R.
, 1999, “
How Heterogeneous Bronchoconstriction Affects Ventilation Distribution in Human Lungs: A Morphometric Model
,”
Ann. Biomed. Eng.
0090-6964,
27
, pp.
14
22
.
25.
Nucci
,
G.
,
Tessarin
,
S.
, and
Cobelli
,
C.
, 2002, “
A Morphometric Model of Lung Mechanics for Time-Domain Analysis of Alveolar Pressures During Mechanical Ventilation
,”
Ann. Biomed. Eng.
0090-6964,
30
, pp.
537
545
.
26.
Olufsen
,
M. S.
, 1999, “
Structured Tree Outflow Condition for Blood Flow in Larger Systemic Arteries
,”
Am. J. Physiol.
0002-9513,
45
, pp.
H257
H268
.
27.
Olufsen
,
M. S.
,
Peskin
,
C. S.
,
Kim
,
W. Y.
,
Pedersen
,
E. M.
,
Nadim
,
A.
, and
Larsen
,
J.
, 2000, “
Numerical Simulation and Experimental Validation of Blood Flow in Arteries With Structured-Tree Outflow Conditions
,”
Ann. Biomed. Eng.
0090-6964,
28
, pp.
1281
1299
.
28.
Steele
,
B. N.
,
Olufsen
,
M. S.
, and
Taylor
,
C. A.
, 2007, “
Fractal Network Model for Simulating Abdominal and Lower Extremity Blood Flow During Resting and Exercise Conditions
,”
Comput. Methods Biomech. Biomed. Eng.
1025-5842,
10
, pp.
39
51
.
29.
Vignon-Clementel
,
I.
,
Figueroa
,
C.
,
Jansen
,
K.
, and
Taylor
,
C.
, 2006, “
Outflow Boundary Conditions for Three-Dimensional Finite Element Modeling of Blood Flow and Pressure in Arteries
,”
Comput. Methods Appl. Mech. Eng.
0045-7825,
195
, pp.
3776
3796
.
30.
Tuminaro
,
M.
,
Heroux
,
S.
,
Hutchinson
,
S.
, and
Shadid
,
J. N.
, 1999, “
Aztec User Guide
,” Version 2.1.
31.
Zamir
,
M.
, 2000,
The Physics of Pulsatile Flow
,
Springer
,
New York
.
32.
Tawhai
,
M. H.
,
Hunter
,
P.
,
Tschirren
,
J.
,
Reinhardt
,
J.
,
McLennan
,
G.
, and
Hoffman
,
E. A.
, 2004, “
CT-Based Geometry Analysis and Finite Element Models of the Human and Ovine Bronchial Tree
,”
J. Appl. Physiol.
8750-7587,
97
, pp.
2310
2321
.
33.
Kitaoka
,
H.
,
Takaki
,
R.
, and
Suki
,
B.
, 1999, “
A Three-Dimensional Model of the Human Airway Tree
,”
J. Appl. Physiol.
8750-7587,
87
, pp.
2207
2217
.
34.
Wall
,
W. A.
,
Wiechert
,
L.
,
Comerford
,
A.
, and
Rausch
,
S.
, 2010, “
Towards a Comprehensive Computational Model for the Respiratory System
,”
International Journal for Numerical Methods in Biomedical Engineering
, in press.
35.
Wiechert
,
L.
, and
Wall
,
W. A.
, 2010, “
A Nested Dynamic Multi-Scale Approach for 3D Problems Accounting for Micro-Scale Multi-Physics
,”
Comput. Methods Appl. Mech. Eng.
0045-7825,
199
, pp.
1342
1351
.
36.
Pedley
,
T. J.
, 1977, “
Pulmonary Fluid Dynamics
,”
Annu. Rev. Fluid Mech.
0066-4189,
9
, pp.
229
274
.
37.
Tawhai
,
M. H.
,
Pullan
,
A. H.
, and
Hunter
,
P. J.
, 2000, “
Generation of an Anatomically Based Three-Dimensional Model of the Conducting Airways
,”
Ann. Biomed. Eng.
0090-6964,
28
, pp.
793
802
.
38.
Shaw
,
R. J.
,
Djukanovic
,
R.
,
Tashkin
,
D. P.
,
Millar
,
A. B.
,
Du Bois
,
R. M.
, and
Corris
,
P. A.
, 2002, “
The Role of Small Airways in Lung Disease
,”
Respir. Med.
0954-6111,
96
, pp.
67
80
.
39.
Evans
,
D. J.
, and
Green
,
M.
, 1998, “
Small Airways: A Time To Revisit?
,”
Thorax
0040-6376,
53
, pp.
629
630
.
40.
Macklem
,
P.
, 1998, “
The Physiology of Small Airways
,”
Am. J. Respir. Crit. Care Med.
1073-449X,
157
, pp.
S181
S183
.
41.
Mols
,
G.
,
Priebe
,
H. J.
, and
Guttmann
,
J.
, 2006, “
Alveolar Recruitment in Acute Lung Injury
,”
Br. J. Anaesth.
0007-0912,
96
, pp.
156
166
.
42.
Kleinstreuer
,
C.
, and
Zhang
,
Z.
, 2009, “
An Adjustable Triple-Bifurcation Unit Model for Air-Particle Flow Simulations in Human Tracheobronchial Airways
,”
ASME J. Biomech. Eng.
0148-0731,
131
, pp.
021007
.
43.
Gemci
,
T.
,
Ponyavinc
,
V.
,
Chena
,
Y.
,
Chen
,
H.
, and
Collins
,
R.
, 2008, “
Computational Model of Airflow in Upper 17 Generations of Human Respiratory Tract
,”
J. Biomech.
0021-9290,
41
, pp.
2047
2054
.
44.
Wiechert
,
L.
,
Metzke
,
R.
, and
Wall
,
W. A.
, 2009, “
Modeling the Mechanical Behavior of Lung Tissue at the Micro-Level
,”
J. Eng. Mech.
0733-9399,
135
, pp.
434
438
.