Abstract
A threedimensional (3D) model of the human airway tree is proposed using a deterministic algorithm that can generate a branching duct system in an organ. The algorithm is based on two principles: 1) the amount of fluid delivery through a branch is proportional to the volume of the region it supplies; and 2) the terminal branches are arranged homogeneously within the organ. These principles define the basic process of branching: generation of the dimensions and directionality of two daughter branches is governed by the properties of the parent branch and the region the parent supplies. The algorithm is composed of nine basic rules and four complementary rules. When the contour of an organ and the position of the trunk are specified, branches are successively generated by the algorithm. Applied to the human lung, the algorithm generates an airway tree that consists of ∼54,000 branches. Its morphometric characteristics are in good agreement with those reported in the literature. The algorithm and the 3D airway model are useful for studying the structurefunction relationship in the lung.
 branching
 threedimensional simulation
 selfsimilarity
 diameter exponent
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 threedimensional (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 (1014, 20,2529). 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
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
In the following, we will assume that each branch is a circular, rigid tube with a constant diameter. An optimum relationship between the flow rate through and diameter of a segment in a living organ was first proposed by Murray (20) and was subsequently used or modified by others (7, 13, 26)
When the flows are conserved before and after branching, we can relate the diameter of the parentd
_{0} to the diameters of the daughtersd
_{1} andd
_{2} (see Fig.1)
On the basis of a minimum energy loss principle, the value ofn was suggested to be three (13, 21). By using various morphometric data, average values ofn 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
Similarly to Eq. 1
, several optimum relationships exist between the branching angles, defined as the angle between the parent and the daughter branch (see Fig. 1), and flow rate in the daughter branches. Using the assumption that the total volume of the three branches in a bifurcation is minimized, Kamiya et al. (13) derived the following equation between the branching angles (θ) and the diameters
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 byEq. 6 .
Relationship Between Length and Diameter
Most morphometric studies have revealed that there is a positive correlation between the diameter and the length of a branch (10, 16,26, 29). By examining the data of Raabe et al. (25), we found that the distribution of the lengthtodiameter ratios showed a Gaussianlike distribution with a mean of 2.8 and an SD of 1.0.
DESIGN OF A BRANCHING DUCTAL SYSTEM
Basic Rules for Generating a Branching Tree Structure
Given a parent branch and the volume region of the parent, we establish the following nine basic rules to generate two daughter branches. Iterative application of these rules to all successive daughter branches results in the generation of a branching duct system within the original volume region of the first parent.
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 “spacedividing plane.” The spacedividing 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 flowdividing ratior is set to be equal to the volumedividing 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 substitutingr 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.
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
Rules 1 and2 reflect the fact that no observations made in real lungs have been reported that would contradict these simple rules.
Rule 3 is based on the assumption that during quasisteady inspiratory flow, the volume increase of the branches is negligible compared with the volume increase of the alveoli.
Rule 4 represents a very simple spacedividing scheme. However, this spacedividing 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 spacedividing 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 onetoone correspondence between a branch and the region it supplies. Therefore, a branch must not leave its own region.
Supplementary Rules
There are four supplementary rules that are introduced to realize a more isotropic space division in the sense defined above. Without the supplementary rules, the shape of a region supplied by a single branch may often become extremely long or flat, which is not observed in real lungs.
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.3 A). On the other hand, if ‖θ_{1} − θ_{2}‖ ≧10°, the division of the parent region into daughter regions, the flowdividing 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 spacedividing halfplanes. The original spacedividing 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 spacedividing halfplane. The second spacedividing halfplane 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. 3 B). These new spacedividing halfplanes define a new volumedividing ratio. Next, basic rules 5 and6 can be reapplied to obtain new diameters and branching angles of the daughter branches. At any given branching, rule 4a is only performed once.
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. ThusEq. 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 ofEq. 6 in this case will not result in the most efficient spacefilling 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 spacedividing plane is determined according torule 4A. We first define “a volumebisecting 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 _{1} andB _{1} for the smaller daughter branch, andA _{2} andB _{2} for the larger daughter branch) by the volumebisecting planes (dashed lines). We also define a volumebisecting angle (ψ in Fig. 4), which is the angle between the parent branch and the volumebisecting plane. Thusrule 6a says that if the angle calculated from Eq. 6 is smaller than the volumebisecting angle, the branching angle is corrected and set to the mean of the volumebisecting angle and the original branching angle. However, if the value of the mean exceeds 90°, the branching angle is set to 90°.
Rule 7a.
The lengthtodiameter ratio is originally set at 3.0. However, it is possible that once the length of a branch is defined by basicrule 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 “distancetolength ratio” as shown in Fig. 5. Although there are no morphometric data on this ratio, many computed tomography images seem to suggest that the distancetolength ratio ranges from 3.0 to 6.0. Thus rule 7aspecifies that the lengthtodiameter ratio, originally set as 3.0, is increased or decreased in steps of 0.25 until the distancetolength ratio falls within its specified lower and upper limits and is not smaller than 1.0.
Rule 8a.
The rotation angle of successive branching planes is originally set at 90° (Fig. 2). However, when the corresponding volumedividing ratio is too small, the daughter branch generated will also be too small. For example, when the volumedividing ratio (and hencer) 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 volumedividing ratio. When the volumedividing 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 volumedividing ratio exceeds this threshold. We regard this threshold of the volumedividing 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 volumedividing 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
Each branch has a set of data including the 3D position of the origination of the branch, the unit vector of the direction of the branch, the flow rate, the normal vector of the branching plane, and sets of planes that define the boundaries of the region supplied by the branch. Additionally, the positions of the ancestor branches are also stored, which is necessary to exclude the space from the region of the branch that has already been occupied by the ancestor branches. The algorithm follows the sequence of the basic rules and then the possible applications of the supplementary rules. If the flow rate of a branch is larger than the flow rate threshold, the branch becomes a parent. When the lengthtodiameter ratio and the rotation angle of branching planes are determined according to rules 7 and 7a andrules 8 and8a, respectively, the data array of two daughter branches is generated. Then, the daughter branches are checked to determine whether their flow rates are larger than the threshold. If the flow rates are smaller than the threshold, the branching stops; otherwise, one or both daughters become a parent. This procedure is then repeated until basic rule 9 is satisfied (the flow rate becomes smaller than the threshold or the branch extends beyond its region). Calculations are done only within a parent region and do not interfere with neighboring branches. However, the procedure is always the same at any generation and independent of the size of the region. Therefore, the branching system so obtained will be a selfsimilar structure (18).
The volume of a daughter region is calculated as follows. Suppose the parent region is given. The x,y, andz 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 spacedividing plane (or two spacedividing planes in case rule 4ais invoked) of the parent region is determined, the volumedividing 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 volumedividing ratio.
Initial Conditions
It is necessary to provide initial conditions specifying the outer boundary of the organ and the configuration of the trunk. To generate models of the airway tree, we examine two types of initial conditions. The first is a volume with a simple surface, which is useful to test the influence of varying parameters on the final structure. The surface is given by the equation
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, lengthtodiameter 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
In the present algorithm, six parameters are specified at the beginning of the calculation: 1) flow rate threshold (Qc), 2) initial value of lengthtodiameter ratio, 3) initial value of the rotation angle of the branching plane,4) lower limit of the distancetolength ratio, 5) upper limit of the distancetolength ratio, and6) the threshold value of the volumedividing ratio. The value of Qc determines the approximate number of TBs (16). We fix the value of Qc to be 0.00006 so as to match the number of termianl branches in the human airway tree. The effects of the other model parameters will be examined separately. The programming was written in C++. Generating a tree with ∼55,000 branches takes ∼90 min on a 200MHz Pentium IIbased personal computer, independent of whether the outer boundary of the organ is simple or complicated.
EVALUATION OF GENERATED AIRWAY TREES
Performance Index (PI)
To find the best combination of the model parameters, it is necessary to evaluate the usefulness and efficiency of the generated tree models as a function of the model parameters. Considering the physiological functions of the airway tree, we first introduce five quantities (efficiency functions) related to certain efficiency properties of the tree and then define a PI to evaluate the overall performance of the algorithm.
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 3cm 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). TheC_{n} is then the ratio of the mean and SD ofn _{t}. This quantity measures the efficiency of the spacefilling capability of a given tree in terms of variations inn _{t}. Because the distribution of TBs within the lung parenchyma is homogeneous,C_{n} should be small. Thus a tree with a smallerC_{n} represents a better design with an increased spacefilling 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 forC _{a} have been estimated to be 0.22 from 130 acini and 0.42 from 209 acini by Kitaoka and Itoh (15) and HaefeliBleuer 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.
The PI is simply the weighted sum of the above five quantities. As we do not know the proper weights of these quantities, all weights are set to be unity. The above quantities or efficiency functions have been defined so that they represent a better structural or functional property of the tree if they take a smaller numerical value. Thus we consider the overall design of a particular tree to be more efficient than that of another tree if the associated PI is smaller.
Effects of Parameter Variations on Tree Design
First, the effects of the supplementary rules are examined, whereas other parameters are kept fixed (ratio of standard length to diameter of 3.0, standard rotation angle of branching plane of 90°, lower and upper limits of distance length ratio of 3.0–6.0, and 0.05 for threshold of volumedividing ratio). The PI, together with the individual efficiency function values for 12 different airway trees (with or without the four supplementary rules), is summarized in Table1. Among these supplementary rules, the most effective is the one that corrects branching angles (rule 6a), and the least effective is the one that corrects the length of a branch (rule 7a). When all supplementary rules are applied, the failure flow rate is decreased from 0.882 (without supplementary rules) to 0.001, and only 14 branches are terminated according the second part of basic rule 9, that is, with flow rate >Q_{c}. The fraction of branches that needed the angle correction of supplementaryrule 6a was 26% of all branches, and the branching angles were increased by <18° in 70% of these branches. Figure 6 shows 2D projection images of the tree structures generated without the supplementary rules (A), with rules 4a, 6a, and7a(B), and with all supplementary rules (C). As Fig. 6 shows, the algorithm without rule 6a is not able to direct branches back toward the root, and the space filling is quite poor. Because the initial conditions are symmetric, the generated tree structures should also be symmetric. However, a slight asymmetry can be seen in Fig. 6 because the volumedividing ratios are calculated by using a grid in the volume that results in slightly asymmetric volume regions with respect to the root of the tree.
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 lengthtodiameter 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 lengthtodiameter 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 Table2. 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 distancetolength 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 Table3, 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 volumedividing 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 volumedividing ratio was increased, PI significantly increased mostly due to increases in Q_{f} andC_{n} and the final tree design was less efficient. Also, notice that there is a positive correlation between the threshold of the volumedividing 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 Table5, 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 lengthtodiameter ratio of 3.0, a rotation angle of branching planes of 90°, a distancetolength ratio between 3.0 and 6.0, and a threshold value of volumedividing 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
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).
DISCUSSION
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 branchingduct 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 selfsimilar 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 exponentn 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 ofn 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 lengthtodiameter 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 airwaybranching 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 NavierStokes 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 volumedividing ratio, but instead flow asymmetry would also be influenced by directionality. Because flow and airway diameter and branching angles are connected via Eqs.46 , 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 flowdividing 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 volumedividing 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
We have introduced a 3D model of the human airway tree down to TBs generated by a deterministic algorithm that incorporates duct branching and space division. The algorithm is relatively simple, yet it yields realistic tree structures with morphometric characteristics that are very similar to those obtained from airway cast studies. The algorithm1) may find useful applications in generating branching structures for computed tomography analysis and for testing various reconstruction algorithms and2) offers a new methodology for a deeper understanding of structurefunction relationships of the tracheobronchial tree such as 3D modeling of airwayparenchymal interactions.
Acknowledgments
We thank Drs. O. G. Raabe and H. C. Yeh and the Lovelace Foundation for providing us with data, which were obtained from research performed at the Inhalation Toxicology Research Institute supported by the National Institute of Environmental Health Science under an Interagency agreement with the US Energy Research and Developmental Administration (now the Dept. of Energy; contract no. DEAC04–76EV01013). We also thank Drs. E. R. Weibel, Bern University; H. Itoh, Kyoto University; K. Eguchi, National Shikoku Cancer Center Hospital; and H. Yamashita for helpful comments and encouragement.
Footnotes

Address for reprint requests and other correspondence: H. Kitaoka, Div. of Functional Diagnostic Imaging, Osaka Univ. Medical School, 2–2 Yamadaoka, Suita City, Osaka 5650871, Japan (Email:kitaokah{at}image.med.osakau.ac.jp).

This study was supported in part by National Science Foundation Grant BES9813599.

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “advertisement” in accordance with 18 U.S.C. §1734 solely to indicate this fact.
 Copyright © 1999 the American Physiological Society