INTRODUCTION

TOP ABSTRACT INTRODUCTION BRANCHING LAWS ESTABLISHED FROM... DESIGN OF A BRANCHING... EVALUATION OF GENERATED AIRWAY... A 3D MODEL OF... DISCUSSION REFERENCES

STRUCTURAL MODELS of the branching ductal systems of an organ are useful for studying the correlation between structure and function. With regard to the lung, Weibel (29) and Horsfield et al. (11) proposed structural airway models that include airway dimensions and connectivity. For this reason then, both of these models have been extensively applied to the study of lung mechanics and gas exchange. However, these models do not include information about the spatial arrangement of the airway structure, and hence they are limited to modeling lung function in one dimension (1D). Because of the rapid advancement of three-dimensional (3D) imaging of the lung that has also gained popularity in the clinical diagnosis of various lung diseases, there appears a strong need for 3D models of the airway and vascular branching systems. For example, the dimensions of central airways have been imaged before and after airway constriction (2). Gillis and Lutchen (5) recently examined the influence of heterogeneous bronchoconstriction on ventilation distribution by using a variation of Horsfield's airway tree model. Had their airway model been arranged in 3D, they would have been able to predict images of ventilation distribution in asthmatic patients with potential clinical impacts. Although there are several geometric airway models such as those proposed by Nelson and Manchester (21) and Martonen et al. (19), these models are still limited to 2D. More recently, Glenny and Robertson (6) and Parker et al. (22) proposed 3D models of the pulmonary arterial system for simulating blood flow distributions in the lung. Although these models have proved useful in predicting function, they are still somewhat limited because, unlike in the lung, the branching structure in these models is symmetric.

Here, we propose a 3D model of branching duct systems that can supply fluid to a given organ evenly distributed in space and apply it to the human tracheobronchial tree. The algorithm that generates the 3D tree structure is based on two simple principles that can be used to determine the branch sizes and spatial arrangement within the organ. The first principle says that the amount of fluid delivery through a branch is proportional to the volume of the region it supplies. This allows us to connect the fluid flow through the branch to the spatial arrangement of the branch. Similar to a lung segment in which the size and the direction of the segmental bronchus are matched to the volume and the shape of the segment, a region in our model and the fluid delivered by the branch supplying that region are matched. This means that the division of the volume of a parent region to daughter regions, and the division of the flow rate in the parent to flow rates in the daughter branches, should be proportional at each bifurcation. This division continues down the tree until the daughter branches become terminal branches. The second principle says that the terminal branches of the tree are homogeneously arranged within the organ. For the human airway tree, the terminal branches of the model are taken to be the terminal bronchioles (TBs) because they are the airway segments beyond which diffusion dominates convective airflow. The pulmonary acini, defined by Fletcher et al. (3) as the regions of the lung distal to the TBs, are then the terminal structural units of the model that uniformly sample the available space and receive nearly identical flows per unit volume. Thus, using our algorithm that is based on the above two principles, effective and homogenous fluid transportation can be realized within an organ.

The organization of the paper is as follows. First, we introduce several rules for the branching geometry that are based on various morphometric studies published in the literature (10-14, 20, 25-29). Next, we describe the basic concepts of the algorithm generating a branching duct system and specify the actual algorithm with conditions appropriate for the human airway tree. We then introduce a method to critically evaluate and compare various special cases of the generated tree structures. Finally, we present the best model of the human airway tree and discuss its limitations and utility for future work.

	BRANCHING LAWS ESTABLISHED FROM PREVIOUS MORPHOMETRIC STUDIES

TOP ABSTRACT INTRODUCTION BRANCHING LAWS ESTABLISHED FROM... DESIGN OF A BRANCHING... EVALUATION OF GENERATED AIRWAY... A 3D MODEL OF... DISCUSSION REFERENCES

From the functional point of view, a simple flow division model can provide efficient fluid transport to the terminal regions or units of a tree structure (16). Therefore, to construct an algorithm for generating an efficient 3D branching system, we need to establish optimal relationships between fluid flow through a segment and the dimensions and directionality of that segment. Thus we first discuss how to connect fluid flow through a segment to the diameter of the segment. Next, we describe two relationships that relate fluid flow division between the daughter segments to the angles of the daughters with respect to their parent. Last, we present arguments for choosing the length of a segment if its diameter and the region it supplies are specified. Once all these relationships are available, we can use them to develop an algorithm that can generate 3D tree structures.

Relationship Between Diameter and Flow Rate

(1)

(2)

When the flows are conserved before and after branching, we can relate the diameter of the parent d₀ to the diameters of the daughters d₁ and d₂ (see Fig. 1)

(3)

We define the flow-dividing ratio (r) as the ratio of flow in the daughter branch receiving the smaller flow to that in the parent. By using Eq. 2, the diameters d₁ and d₂ of the two daughter branches can now be expressed in terms of d₀ and r as

(4)

where we require that 0 < r  0.5 and d₁  d₂.

View larger version (9K): [in this window] [in a new window]   Fig. 1.   Dichotomous branching that occurs in a single plane. Diameters of parent and 2 daughters are d0, d1, and d2, respectively. Flow-dividing ratio r is ratio of flow rate in daughter branch receiving smaller flow to flow rate in parent. Branching angles 1 and 2 are related to r via Eq. 6.

On the basis of a minimum energy loss principle, the value of n was suggested to be three (13, 21). By using various morphometric data, average values of n have been reported by Groat (7) to be 2.6 for arteries, by Suwa et al. (26) to be 2.7 for arteries, and by Horsfield and Thurlbeck (12) to be 2.4-2.9 for airways of four mammals. The reason for this discrepancy between the theoretical value and the morphometric values is still unclear. Recently, we examined and rederived Eq. 3 by using morphometric data of four mammalian airway trees published by Raabe et al. (25), and we estimated n to be 2.8 in human lungs (16).

Relationship Between Branching Angle and Flow Rate

(5)

(6)

Measurement of branching angles is difficult so that there is a paucity of angle data in the literature. Moreover, there are no data comparing angles and flow rates in airway or vessel segments. Nevertheless, we can assess the approximate validity of Eq. 5 by examining some morphometric data. Thurlbeck and Horsfield (27) investigated the dependence of branching angles on the difference between branching orders at a bifurcation, which should reflect the difference between flow rates carried by the two daughter branches. They found a tendency similar to that predicted by Eq. 6.

Relationship Between Length and Diameter

	DESIGN OF A BRANCHING DUCTAL SYSTEM

TOP ABSTRACT INTRODUCTION BRANCHING LAWS ESTABLISHED FROM... DESIGN OF A BRANCHING... EVALUATION OF GENERATED AIRWAY... A 3D MODEL OF... DISCUSSION REFERENCES

Basic Rules for Generating a Branching Tree Structure

Rule 1: Branching is dichotomous.

Rule 2: The parent branch and its two daughter branches lie in the same plane, called the branching plane.

Rule 3: The volumetric flow rate through the parent branch is conserved after branching; that is, the sum of the flows in the daughter branches is equal to the flow in the parent branch.

Rule 4: The region supplied by a parent branch is divided into two daughter regions by a plane called the "space-dividing plane." The space-dividing plane is perpendicular to the branching plane and extends out to the border of the parent region. There is a supplementary rule (rule 4a) for correcting the space division scheme whenever the shape of the region requires it.

Rule 5: The flow-dividing ratio r is set to be equal to the volume-dividing ratio, defined as the ratio of the volume of the smaller daughter region to that of its parent.

Rule 6: Diameters and branching angles of the two daughter branches are determined by substituting r from rule 5 into Eqs. 4 and 6, respectively. There is a supplementary rule (rule 6a) for correcting the branching angle according to the shape of the daughter region.

Rule 7: The length of each daughter branch is assigned a value that is three times its diameter. There is a supplementary rule (rule 7a) for correcting the length according to the shape of the region.

Rule 8: If branching continues in a given direction, the daughter branch becomes the new parent branch, and the associated branching plane is set perpendicular to the branching plane of the old parent (Fig. 2). We call the angle between the two successive branching planes the "rotation angle of the branching plane." There is a supplementary rule (rule 8a) for correcting the rotation angle according to the shape of the region.

View larger version (16K): [in this window] [in a new window]   Fig. 2.   Rotation angle of successive branching planes. Branching plane of daughter branch is initially set to be perpendicular to branching plane of parent.

Rule 9: The branching process in a given direction stops whenever the flow rate becomes less than a specified threshold or the branch extends beyond its own region.

Explanations of the Basic Rules

Rule 3 is based on the assumption that during quasi-steady inspiratory flow, the volume increase of the branches is negligible compared with the volume increase of the alveoli.

Rule 4 represents a very simple space-dividing scheme. However, this space-dividing scheme is corrected according to supplementary rule 4a, if necessary (see rule 4a below).

Rule 5 is a natural consequence of our hypothesis that the main function of the branching structure is to supply similar amounts of fluid to the acini that uniformly fill in the space of the organ.

Rule 6 expresses the optimal relationships between flow rate and diameter and flow rate and branching angles, as specified by Eqs. 4 and 6, respectively.

Rule 7 simply determines the branch length in accordance with our observations on the relationship between airway length and diameter from the data of Raabe et al. (25).

Rule 8 sets up the conditions for the next branching. There have been no reports on the rotation angle of branching planes. However, close examinations of airway casts and bronchographic images suggest that the rotation angle of the branching plane is close to a right angle. According to rule 4, space division occurs perpendicularly to the branching plane. The choice of 90° for the rotation angle of the space-dividing planes is advantageous if the initial shape of the volume (to be filled in by the tree structure) is not very irregular. In this case, successive application of perpendicular branching planes leads to a reasonably isotropic space division. By isotropic space division we mean a subdivision of a given volume into many smaller volume elements that have similar shape and volume in all directions.

Rule 9 reflects two conditions that stop branching. The first says that fluid transport changes from conductive flow to diffusion in the airways or capillary network transport in the vasculature when branching reaches the terminal segments of the tree. The simplest realization of this concept is to set up a flow threshold value. The second condition says that it is necessary to assume a one-to-one correspondence between a branch and the region it supplies. Therefore, a branch must not leave its own region.

Supplementary Rules

Rule 4a. The purpose of this supplementary rule is to create a more isotropic space division than is possible with basic rule 4. Once the two branching angles are calculated by using basic rule 6, they are examined for the following conditions. If the two branching angles satisfy the relationship |_1  2| < 10°, there is no change in the angles (Fig. 3A). On the other hand, if |_1  2| 10°, the division of the parent region into daughter regions, the flow-dividing ratio r and, consequently, the diameters and angles of both daughters are not optimal and are reevaluated as follows. First, the parent region will be divided into two new daughter regions with two separate space-dividing half-planes. The original space-dividing plane that included the parent branch and extended out to the boundary is now cut at the tip of the parent branch. This plane defines the first space-dividing half-plane. The second space-dividing half-plane is also perpendicular to the branching plane and starts at the tip of the parent and extends out to the border of the available space, in the direction of the bisector angle of the two daughter branches (see Fig. 3B). These new space-dividing half-planes define a new volume-dividing ratio. Next, basic rules 5 and 6 can be reapplied to obtain new diameters and branching angles of the daughter branches. At any given branching, rule 4a is only performed once.

View larger version (18K): [in this window] [in a new window]   Fig. 3.   Branching plane showing parent and 2 daughter branches and how space-dividing planes divide parent region into 2 daughter regions. A: when difference between branching angles is <10°, a single plane (thin line) divides parent region into daughter regions. B: when difference between branching angles is 10°, space division is done by 2 separate space-dividing planes: 1 that contains parent and is perpendicular to branching plane and 1 that connects end tip of parent and lower right boundary of parent region (thin line). These new space-dividing half-planes define a new volume-dividing ratio. Basic rules 5 and 6 are reapplied to obtain new branching angles of daughter branches ('1 and '2).

Rule 6a. According to Eq. 6, the maximum branching angle is 90°. When the smaller daughter branch serves a small region, r will take a small value. For example, for r = 0.01, the smaller angle is 78 and 80° for n = 2.8 and n = 3, respectively. Thus Eq. 6 predicts that practically all branching angles should be <80°. However, it is generally accepted that there are branches with a branching angle that is a nearly right angle, called a monopody (29). This situation occurs when the region supplied by a branch extends backward in the direction of the parent branch, as in the superior segmental branches of bilateral lower lobes in the human lung. It is obvious that application of Eq. 6 in this case will not result in the most efficient space-filling arrangement of the descendent branches within the region of the parent because, in reality, those branching angles are larger than the values predicted by Eq. 6. Thus the purpose of rule 6a is to correct branching angles in these situations. It is applied after the space-dividing plane is determined according to rule 4A. We first define "a volume-bisecting plane," which includes the branching point, is perpendicular to the branching plane, and divides the daughter region into two regions with equal volumes. For example, both daughter regions in Fig. 4 are divided into two equal volumes (A₁ and B₁ for the smaller daughter branch, and A₂ and B₂ for the larger daughter branch) by the volume-bisecting planes (dashed lines). We also define a volume-bisecting angle ( in Fig. 4), which is the angle between the parent branch and the volume-bisecting plane. Thus rule 6a says that if the angle calculated from Eq. 6 is smaller than the volume-bisecting angle, the branching angle is corrected and set to the mean of the volume-bisecting angle and the original branching angle. However, if the value of the mean exceeds 90°, the branching angle is set to 90°.

View larger version (20K): [in this window] [in a new window]   Fig. 4.   Correction of branching angles to achieve space-filling arrangement according to supplementary rule 6a. Volume-bisecting planes (dotted lines) include branching point, are perpendicular to branching plane, and divide daughter regions into 2 regions with equal volumes (A1 = B1 in region of smaller daughter branch and A2 = B2 in region of larger daughter branch). Volume-bisecting angle 1 is angle between parent branch and volume-bisecting plane. Supplementary rule 6a says that if angle calculated from Eq. 6 is smaller than volume-bisecting angle, branching angle is set to mean of volume-bisecting and original branching angles.

Rule 7a. The length-to-diameter ratio is originally set at 3.0. However, it is possible that once the length of a branch is defined by basic rule 7, the branch will extend beyond the parent region or its end will come too close to the boundary of the region. To avoid this situation, we define a "distance-to-length ratio" as shown in Fig. 5. Although there are no morphometric data on this ratio, many computed tomography images seem to suggest that the distance-to-length ratio ranges from 3.0 to 6.0. Thus rule 7a specifies that the length-to-diameter ratio, originally set as 3.0, is increased or decreased in steps of 0.25 until the distance-to-length ratio falls within its specified lower and upper limits and is not smaller than 1.0. 

View larger version (10K): [in this window] [in a new window]   Fig. 5.   Definition of distance-to-length ratio on branching plane. P0, starting point of the branch; P1, ending point of the branch; P2, cross point between a prolongation of the branch and the boundary of the region.

Rule 8a. The rotation angle of successive branching planes is originally set at 90° (Fig. 2). However, when the corresponding volume-dividing ratio is too small, the daughter branch generated will also be too small. For example, when the volume-dividing ratio (and hence r) is 0.05, the diameter of the daughter branch as obtained from Eq. 4 is only ~30% of the parent diameter. Changes in diameter that are so large are rarely observed in real lungs. Therefore, we set up a threshold value for the lower limit of the volume-dividing ratio. When the volume-dividing ratio for a rotation angle of 90° is less than this threshold, the rotation angle of the branching plane is increased or decreased in steps of 9° until the volume-dividing ratio exceeds this threshold. We regard this threshold of the volume-dividing ratio as a parameter in the model with a value of 0.05. Additionally, when the flow rate in a branch becomes very small, <1.5 times the flow rate threshold, the lower limit of the volume-dividing ratio is increased to 0.35. This will guarantee that the terminal volumes (e.g., acini for the airway tree) do not become unrealistically small.

Data Structure and Algorithm

The volume of a daughter region is calculated as follows. Suppose the parent region is given. The x, y, and z coordinates of the parent region occupied in the 3D space are determined. Next, a rectangular space including the entire parent region is set up and divided into small rectangles with sides that are 20 times smaller than those of the embedding rectangle. When the space-dividing plane (or two space-dividing planes in case rule 4a is invoked) of the parent region is determined, the volume-dividing ratio is calculated as a ratio of the number of rectangles in the daughter region to that in the parent region. The volume of the daughter region is simply the product of the volume of the parent region and the volume-dividing ratio.

Initial Conditions

The second set of initial conditions specifies a realistic boundary mimicking the shape of the lung and is given by a combination of the following approximate surfaces: the thorax modeled as a parabolic ellipsoid surface, the thoracic vertebrae as a sinusoidal surface, the diaphragm as a parabolic ellipsoid surface, the heart as an ellipsoid, the descending aorta as a cylinder, and all other mediastinal parts as simple planes. Branching angles, length-to-diameter ratio, and rotation angle of each branching plane down to segmental branches are assigned in advance so as to meet the standard branching pattern of the central airways reported by Yamashita (31). Maximum width of the inner surface of the thorax is taken to be 30 cm, and the diameter of the trachea is 18 mm. Total lung volume is then estimated to be 6,388 ml, consistent with that in normal adults.

Model Parameters

	EVALUATION OF GENERATED AIRWAY TREES

TOP ABSTRACT INTRODUCTION BRANCHING LAWS ESTABLISHED FROM... DESIGN OF A BRANCHING... EVALUATION OF GENERATED AIRWAY... A 3D MODEL OF... DISCUSSION REFERENCES

Performance Index (PI)

Failure flow rate (Q_f). The terminal branches of a tree that are terminated because they extend beyond their own region (second part of basic rule 9) violate the concept that fluid should be delivered homogeneously in the space available for the tree. These branches are called "failed terminal branches," and Q_f is calculated as the sum of the flow rates carried by all failed terminal branches. Because the total flow rate is normalized to unity, 100% × Q_f is the percent error in optimal fluid delivery.

Coefficient of variation of the density of terminals (C_n). The lung volume is divided into 181 cubes, each having a 3-cm side. For our lung shape defined above, 75% of these cubes are completely within the contour of the lung. The total number of terminal branches (NT) generated by the algorithm is counted in each of the 181 cubes. NT is then normalized by the volume of the cubes (i.e., 27 ml), which provides the density of terminals, n_t, in units of terminal number per unit volume (in milliliters). The C_n is then the ratio of the mean and SD of n_t. This quantity measures the efficiency of the space-filling capability of a given tree in terms of variations in n_t. Because the distribution of TBs within the lung parenchyma is homogeneous, C_n should be small. Thus a tree with a smaller C_n represents a better design with an increased space-filling capability.

Coefficient of variation of acinar volumes (C_a). We define the acinus here as the terminal region supplied by a terminal branch (i.e., with flow rate <Qc). C_a is defined as the coefficient of variation of the volumes of the terminal regions. The value of C_a should be small for uniform space division. Values for C_a have been estimated to be 0.22 from 130 acini and 0.42 from 209 acini by Kitaoka and Itoh (15) and Haefeli-Bleuer and Weibel (8), respectively.

Volume ratio (V_r). V_r is defined as the ratio of the total airway volume of the model to the volume of the space within the lung contour. V_r is simply the dead space of the airway tree relative to lung volume and should take small values (<0.1).

Energy loss per terminal (EL). The energy loss associated with fluid flowing through a single branch is defined as the fluid flow resistance of the branch multiplied with the square of the flow rate. The total energy loss for delivery of fluid to the acini is the sum of the individual energy losses over all branches. Because the total energy loss increases with the number of terminals, we compare trees using the total energy loss normalized with the number of terminals. Although the previous quantities assessed the efficiency of the structure of a tree, EL is a measure of the function of the tree, i.e., the efficiency with which fluid can be delivered to the acini. In the actual calculations of EL, we assume Poiseuille flow and, for simplicity, the diameter of the root and the viscosity of the fluid are taken to be unity.

Effects of Parameter Variations on Tree Design

View larger version (37K): [in this window] [in a new window]   Fig. 6.   Projected images of tree structures constructed with simple lung contour. A: without supplementary rules. B: with rules 4a, 7a, and 8a. C: with all supplementary rules.

The above results demonstrate that the supplementary rules substantially improved the design of the tree. Thus, in an examination of the effects of other parameters, all supplementary rules were kept in the design. The mean and SD of the length-to-diameter ratios that the algorithm realized by the application of rule 7a were 3.06 ± 0.46, close to the values obtained from data of Raabe et al. (25) as mentioned above. When the initial value of the length-to-diameter ratio was changed to 2.7 or 3.3, the values realized by the algorithm and the PI did not change appreciably.

The mean and SD of the rotation angle of the branching plane realized by the algorithm (using supplementary rule 8a) were 90 ± 7°. Changing the initial value of the rotation angle from 90° led to a significant increase in Q_f and PI, as shown in Table 2. Thus the choice of a right angle for the rotation angle of the branching plane appears to be advantageous.

Next, the effects of the lower and upper limits of the distance-to-length ratio were examined by varying the lower limit from 2.5 to 4.5 and the upper limit from 4.0 to 8.0. As shown in Table 3, the combination of 3.0 and 6.0 was the best, although the effect of these parameters on the overall tree design remains small.

The effect of the threshold value of the volume-dividing ratio is summarized in Table 4. When the threshold was lowered to 0.02, Q_f decreased but at the expense of an increased coefficient of variation of acinar volumes. When the volume-dividing ratio was increased, PI significantly increased mostly due to increases in Q_f and C_n and the final tree design was less efficient. Also, notice that there is a positive correlation between the threshold of the volume-dividing ratio, the SD of realized rotation angle, and Q_f.

Although the diameter exponent was not regarded as a model parameter, we examined its effect on the tree design. As shown in Table 5, this parameter can have a significant effect on many features of the tree, including total airway volume and energy loss. Specifically, for n between 2.8 and 3.0, there are only slight differences in the spatial arrangement of the tree structures, but there is a significant difference in properties related to function, i.e., energy loss. This issue warrants further consideration.

In summary, the best conditions for designing a tree in a shape similar to the lung are given by the combination of the following values of the parameters: a length-to-diameter ratio of 3.0, a rotation angle of branching planes of 90°, a distance-to-length ratio between 3.0 and 6.0, and a threshold value of volume-dividing ratio of 0.05. Finally, we also examined effects of the boundary conditions on tree design by using other geometric shapes such as a sphere, cylinder, and cube. The results were similar to those obtained with our lung shape. This suggests that the algorithm is robust and capable of generating an efficient tree structure into many different organ shapes.

	A 3D MODEL OF THE HUMAN AIRWAY TREE

TOP ABSTRACT INTRODUCTION BRANCHING LAWS ESTABLISHED FROM... DESIGN OF A BRANCHING... EVALUATION OF GENERATED AIRWAY... A 3D MODEL OF... DISCUSSION REFERENCES

On the basis of the above results, we constructed a 3D model of the human airway tree by using all supplementary rules and the second, more complicated boundary shape. The corresponding five efficiency function values were as follows: Q_f = 0.018, C_n = 0.18, C_a = 0.34, V_r = 0.027, EL = 0.011, and PI = 0.576. Figure 7 shows an anterior and a right lateral view of this 3D model. Different branches in Fig. 7 distal to different segmental bronchi are shown by different colors. The tree structure itself is remarkably similar to the real airway tree. There are 27,306 terminal branches in the model that correspond to TBs. The mean and SD of the generation numbers of TBs are 17.6 ± 3.4, with a minimum value of 8.0 and a maximum value of 32.0. The distribution of the generation numbers of the TBs is wider and the mean value is higher than those for the model of Horsfield et al. (11). The reason is most likely due to the more frequent asymmetric branching in our model. If the tracheal diameter is 18 mm, the mean diameter of TBs in the present model is 0.48 ± 0.06 mm, which is slightly smaller than previously reported values (10, 17, 30). The reason, again, seems to be due to the more asymmetric branching pattern in the present model compared with previous models. The lung volume of the model is 6,388 ml, and the corresponding total airway volume is estimated to be 174 ml. This value is consistent with that in Weibel's model (29). The lung parenchyma is ~90% of the total lung volume according to Weibel. The mean volume of the acini is then estimated to be 0.211 ± 0.076 ml, which is consistent with previously reported morphometric data (4, 9, 15, 24).

View larger version (79K): [in this window] [in a new window]   Fig. 7.   Image of 3-dimensional airway tree model obtained with best choice of parameters. It has 54,611 branches, where branches distal to different segmental bronchi are shown in same color as segmental bronchus. A: anterior view. B: right lateral view.

	DISCUSSION

TOP ABSTRACT INTRODUCTION BRANCHING LAWS ESTABLISHED FROM... DESIGN OF A BRANCHING... EVALUATION OF GENERATED AIRWAY... A 3D MODEL OF... DISCUSSION REFERENCES

The purpose of this paper was to develop an algorithm capable of generating a branching tree structure in an arbitrary shape. Applying it to the lung resulted in a 3D model of the human airway tree that is virtually indistinguishable from the real airway tree. The morphometric features of the model are also in good agreement with those obtained from airway casts, suggesting that the algorithm proposed here is effective in growing a branching-duct system in the lung. The algorithm is also applicable to other duct systems, provided the system is large enough to contain a sufficient number of generations. The algorithm is based on a collection of rules that are iteratively applied to successive branches. Because the actual computations do not depend on the size of the region, except for branches that are very close to the boundaries of the organ, the tree exhibits self-similar properties. The 3D tree generated by our algorithm is therefore a finite fractal structure (14, 17, 18).

As we have shown earlier, the distribution of diameters is closely related to the distribution of flows through the segments, both being a power law (16). The exponent in the cumulative power law distribution of diameters is equal to the diameter exponent n in Eq. 1. According to Uylings (28), n depends on whether the flow is laminar or turbulent. The flow in the central airways can be turbulent or transitional, implying a smaller n (<2.7), whereas it is certainly laminar in the periphery of the lung so that n is expected to be close to 3.0. However, in this study we used a single value of n for the sake of simplicity of the algorithm. Further studies both theoretical and experimental are needed to examine this issue in more detail.

The computational process in our algorithm is completely deterministic. The advantage of a deterministic algorithm is that it can generate tree structures reproducibly. One may even be tempted to associate certain steps of the algorithm with various biochemical and molecular signals that initiate the bifurcation of a branch (23). The disadvantage of the deterministic nature of the algorithm is that, in its present form, it is not able to introduce randomness into the structure. Thus the algorithm would not be able explain interindividual variations in airway tree structures within a given species. However, this can be easily remedied by allowing some randomness in certain parameters. When a small amount of fluctuation was added to the initial value of the length-to-diameter ratio (SD = 0.3), the generated trees were slightly different, but the statistical properties of all features of the trees were identical. Also, with the addition of a small number of fluctuations (SD = 5°) to the rotation angle of the branching planes of 90°, the mean of the realized rotation angles and the PI did not change significantly. Another possibility is to add some randomness to the initial conditions when trees are repeatedly generated for a given species. This perturbation also leads to different structures because the branching process and the generated tree structure are quite sensitive to the initial conditions. For example, the direction and length of the trunk or small alterations in the initial shape can have a significant impact on the subtended structure because the small differences in the structure are amplified at successive branching steps. However, the algorithm has shown a remarkable tolerance against these perturbations in the sense that the statistical properties (e.g., efficiency functions, distribution of diameter, and so on) of the tree structures generated under slightly different conditions were the almost identical. Thus small variations in shape and trunk properties or fluctuations in certain parameters can lead to structures that can mimic interindividual differences within a species.

It is well known that different species have different airway-branching patterns. For example, dogs have an airway tree that is very asymmetric compared with a human airway tree, even though the initial size and shape of the thoracic cavity of dogs and humans are not grossly different. Our algorithm is not able to account for such an interspecies difference in structure, and hence additional factors may have to be taken into account. One possibility is that fluid flow inertia can make flow distribution very inhomogeneous and direction dependent. Recently, Andrade et al. (1) solved the Navier-Stokes equation in a 2D airway system and found that the flow distribution coming out of a symmetric tree largely depended on the Reynolds number. For higher Reynolds numbers, convective acceleration produced a highly asymmetric flow pattern where, due to fluid flow inertia, segments that were more aligned with their grandparents received more flows. Thus flows in the central airways would not always follow the volume-dividing ratio, but instead flow asymmetry would also be influenced by directionality. Because flow and airway diameter and branching angles are connected via Eqs. 4-6, the resulting branching structure would also depend on the average Reynolds number of a species. Because the breathing rate in dogs is higher than in humans whereas airway diameters are similar, the average Reynolds number can be higher in dogs. As a consequence, in dogs, flow asymmetry would be more pronounced and hence the corresponding structure would also be more asymmetric. Thus, to account for interspecies differences, specifying different initial conditions may not be sufficient. Indeed, in our previous 1D airway model (16), it was necessary to directly implement different mean flow-dividing ratios to generate human and dog airway tree models.

The present model also has some shortcomings. As the contour of the lung is modeled more and more faithfully, the model provides a larger and larger Q_f. The reason is likely to be the existence of concave surfaces in the real lung contour due to the heart and the aorta. In this case, branching may stop too early according to rule 9 even if flow rates are much larger than Q_c.

As mentioned earlier, the branching pattern in our model is more asymmetric than in the model of Horsfield et al. (11). If a more symmetric branching pattern is needed, the threshold value of the volume-dividing ratio should be increased. Unfortunately, this would result in a significantly increased Q_f (see Table 4). To compensate for this effect, the algorithm would have to be refined. Another shortcoming of the model is related to the branches that are located near and reach back toward one of their ancestors. The number of such branches is higher in the model than in the real lung. In reality, a region neighboring a proximal branch may be supplied by branches other than its descendants. For example, the region may be supplied by a branch that has a shorter total pathway from the trachea. Thus ventilation distribution in the present model may show an increased heterogeneity compared with measured ventilation distribution due primarily to these longer pathway lengths. To solve this problem, it would be necessary to change the rules for space division. However, this could be done only at the expense of a much more complicated algorithm. Nevertheless, the fraction of these branches in the whole tree is small, so that their influence on function may not be important.

Conclusion