An Intelligent Parallel Method for Automatic Production of the Human Bronchial Tree Fullling Bronchopulmonary Segmental Isolation Based on Stochastic Parametric L-system

10 In the human lung, air is transported from the environment to the respiratory zone by conducting 11 airways containing ~60000 airways with different scales. An accurate model of this complicated 12 structure is crucial for studying transport phenomena in the human lung. In the present study, the 13 parametric Lindenmayer System (L-System) has been used to produce the human bronchial tree 14 in the bronchopulmonary segments whose shapes have been obtained based on the medical image 15 processing techniques. The airways propagate into the host segments simultaneously similar to 16 what happens in real life using multi-threads parallelism method. Each module in rewriting rules 17 of the formulated L-System plays the role of an airway with its characteristic properties that have 18 been found almost independently based on an intelligent procedure. This feature enables the 19 method to be less sensitive to change in the host shape and size of the lungs. In addition, the 20 stochastic behavior of the proposed procedure generates slightly different structures with the same 21 statistical properties. Furthermore, the dimensions of the terminal branches are functions of age, 22 which make the proposed method capable of generating the tree structure for the subjects with 23 various ages fulfilling almost all anatomical and physiological requirements. The morphometric 24 characteristics of the generated structure are in good agreement with the corresponding 25 experimental data reported in the literature. The results show that the proposed method 26 outperforms the previously reported models. The accurate three-dimensional model proposed by 27 this study can be used for simulation of fluid distribution, mass and heat transfer, drug delivery, 28 and particle deposition in the human respiratory system. 29


Introduction 31
There are three and two lobes in the right and left lungs, respectively. Each lobe is also divided 32 into sections called bronchopulmonary segments. Layers of connective tissue separate the 33 segments from each other. Each segment is supplied by a segmental bronchus and bronchial artery 34 and functions in a manner that does not depend on the other segments [1]. The right lung consists 35 of 10 segments whereas the number of left lung segments varies from 8 to 10. The most common 36 structure of the left lung comprises 8 segments. Air goes into trachea from nose or mouth during 37 breathing. Then it passes through airways into the lungs and backs out again. The airways in the 38 human lungs can be classified into two main parts: conducting airways and respiratory airways. 39 Conducting airways start from the trachea and propagate into the lung as terminal bronchioles. 40 They supply air from the trachea to respiratory units where gas exchange between respiratory 41 airways and blood takes place. The airways have been confined in a limited space resulting in a 42 complicated structure with multi-scales. 43 The conducting airways of the human lungs form a complex asymmetric tree. Mathematical 44 models of the tracheobronchial tree have been widely used to study aerosol deposition [2], gas 45 transport [3][4][5][6], ventilation distribution [7,8], and heat transfer [9]. The accurate modeling of all 46 aforementioned subjects requires appropriate and reliable modeling of the bronchial tree. For this 47 reason, modeling of the bronchial structure has been attracted much attention. Since the problem 48 does not have a unique answer, the previous research tried to present mathematical models with 49 various simplifications in such a way that the statistical properties of the final structure satisfy the 50 anatomical data. Weibel[10] and Yeh and Schum[11] respectively proposed one and five typical 51 symmetric paths' models of the airways obtained based on cast measurement. Although this 52 approach has the lowest computational cost for studying the aforementioned applications it is not 53 appropriate for studying the effect of anatomical asymmetry. 54 To study the properties of the bronchial tree, the branches should be classified. There are three 55 methods of classification that are commonly used: generation, Horsfield and Strahler order. The 56 first one was used by Weibel[10]. In this approach, counting starts at the trachea, considered as 57 generation 0 or 1, and increases by one in each bifurcation. Strahler and Horsfield orders [12] are 58 the most common classifying methods for asymmetric trees in which the branches are numbered 59 upward from end branches instead of downward starting with the trachea. Horsfield et.al[13] used 60 parameter delta, the difference between the order of the daughters' branches, to express the 61 asymmetry at dichotomy. This method assumes a uniform degree of asymmetry throughout the 62 tree but the branches with the same order have different diameters or lengths in reality. Koblinger 63 and Hofmann [14] proposed a stochastic asymmetric model for human bronchial tree based on 64 morphometric data of Raabe et.al. [15]. However, all mentioned simplified geometries do not take 65 into account the spatial positions of airways. 66 Kitaoka et.al[16] used nine basic rules and four complementary rules to produce a three-67 dimensional tree (3-D) representing the bronchial tree of the lungs. They used a power-law 68 relationship between the flow dividing ratio and diameter. To find the branching angle, the relation 69 proposed by Murray[17] was employed, which has been derived based on the minimization of 70 power dissipation and is a function of the flow dividing ratio. They assigned a value to the branch 71 length that was three times its diameter. Their deterministic model was very sensitive to the initial 72 value of model parameters and unable to generate various structures with identical structural 73 properties. Tahwai et.al [18] proposed a model for the construction of a 3-D bronchial tree in each 74 lobe based on Monte Carlo method. They distributed seed points on the host volume and space 75 boundaries of the bronchopulmonary segments are only manually found once. For other subjects, 145 they can be estimated using the lobes' boundaries assuming that the relative position of the 146 segments' boundaries in the corresponding lobes remain unchanged. 147

Mathematical formulation 148
Since the problem belongs to the category of inverse problems and there is no one to one 149 mapping between the model parameters and the observations, the problem has no unique solution. 150 In what follows, we present various models and show that the model with the least number of 151 parameters is suitable because there is not enough data to determine many parameters. 152 The human bronchial tree can be described by a simple directed graph1 whose edges are three-153 dimensional. All nodes and edges of the graph should be determined in such a way that the 154 structure satisfies all anatomical and physiological constraints of a living organ. Therefore, the 155 problem can be defined by: 156 (2.1) ( ( , )): × → ℝ + :̃( ( , )) <̃ where ( , ) is a simple directed graph with a set of nodes and a set of three dimensional 157 directed edge E. and ̃ are the norms of the vector of objective functions ̃, and vector of 158 constraints whose elements are members of ℝ + . 159

Growth based formulation 160
There are about 60000 conducting airways in the human bronchial tree [23] and thus the number 161 of decision variables of the problem is 240000 calculated by considering four decision variables 162 (i.e. the diameter and the endpoint position) for each branch. It is difficult to find such a large 163 number of variables in such a way that the proposed algorithm remains robust in the presence of 164 small disturbances or deformations. Therefore, it is essential to reduce the number of decision 165 variables in such a way that the model remains accurate and predictable. For this purpose, the 166 graph G can be obtained using a rule-based method such as Lindenmayer system. This method and 167 its application in the modeling of the structure are described in the rest of this section. 168 Lindenmayer system or L-system was introduced and developed by Aristid Lindenmayer[24,169 25]. It is a type of formal grammar and a parallel rewriting system that uses production rules or a 170 set of grammars to produce branching structures. L-systems have been used to describe the growth 171 process of plants [26,27], model morphology of organisms [22,28], and the production of self-172 similar fractals. It can be defined as a triplet L( , , ). L-System consists of four main parts: 173 • : A set of symbols that are used to make string. It consists of both variables, elements that 174 can be replaced, and terminal or constants, the elements that cannot be replaced. The set of 175 all words and set of all nonempty words are denoted by * and + , respectively. 176 • + : An axiom that is an initial string from which the branching starts. 177 • ⊂ × * : A set of production rules that are employed to make strings from each 178 symbol. It consists of a predecessor and a successor and is written in the following format. 179 The symbol "→" is used to separate the predecessor and successor. 180 • Turtle interpretation: a mechanism that interprets the strings into geometric shapes. 181 Lindenmayer used brackets "[" and "]" to produce axial trees. These brackets are employed to 182 enable L-system to produce branching structure. The brackets "[" and "]" are not included in the 183 set V, therefore the set V is extended to set = ⋃{[, ]} for bracketed L-system. 184 If the graph representing the bronchial tree is obtained based on L-system, the formulation of the 185 problem will be: 186 To generate the structure of the human bronchial tree using the above approach, the number of 187 production rules, i.e. decision variables, is numerous because the bronchial tree has a non-uniform 188 structure. Therefore, it is of interest to reduce the number of production rules. For this purpose, the 189 parametric L-system is used to generate the tree with a minimum number of decision variables that 190 include both production rules and parameters. 191 In parametric L-system predecessors and successors are parametric strings with letters from set 192 V and formal parameters, real-valued parameters appearing in words, from set Σ. Any module 193 ( 1 , 2 , … , ) with letter and parameters 1 , 2 , … , ℝ belongs to set × ℝ * . 194 Parametric L-system is denoted by quadruplet L( , Σ, , ) where: 195 197 where ( × ℝ * ) + and ℂ(Σ) are a set of nonempty modules and a set of logical expressions, 198 respectively. ℰ(Σ) denotes the set of all procedures for finding parameters. Symbol "*" above any 199 set B denotes the set of all finite combinations of members of B. The parametric production rule 200 is written in the following format. 201 Where the symbol ":" is used to separate the condition and predecessor. Any formal parameters 202 may appear in the conditions. The reader is referred to the reference [24] for further details on 203 different types of L-System. Using parametric L-system, Eq. (2.3) can be reduced to: 204 (2.5) (L( , Σ, , )) : ̃( L( , Σ, , )) <̃ All structures produced by the same deterministic L-system are identical. To generate the bronchial 205 structure for various people, it is necessary to randomize parameters so that it considers case-to-206 case variations and preserves statistical properties of the structure. In the next section, the process 207 of generating tracheobronchial structure using parametric L-system with random parameters is 208 presented. 209

Application of Parametric L-System in the Modelling od Conducting Airways 210
Trachea branches into main bronchi which are then branching into lobar (secondary) bronchi. 211 Lobar bronchi further divide into segmental (tertiary) bronchi which branch into smaller airways 212 called bronchiole connected to terminal bronchiole feeding respiratory zone. Conducting airways 213 can be divided into two main sections. The first section extends from the trachea to segmental 214 bronchi. The second section comprises segmental bronchi to terminal bronchioles. All branches 215 with identical segmental bronchus are generated into their corresponding bronchopulmonary 216 segment and the generation of branches in each segment occurs in parallel. Parametric L-system 217 is used to produce the branches of each mentioned section. Each airway can be determined using 218 its parent, diameter, length, and orientation. Parameters of parametric L-system should be 219 determined by an intelligent procedure to satisfy anatomic constraints for all cases. The decision 220 variables of the proposed procedure should be determined in such a way that the objective function 221 of Eq. (2.5) is minimized. Each member of the vector ̃ corresponds to the difference between the 222 statistical property of simulated structure and its reference value according to the following Eq. 223 (3.1). 224 Where, , is a function that calculates the corresponding property of the simulated structure. 225 The list of geometrical properties considered in Eq. (3.1) is given in

Central airways 238
As mentioned before, the number of bronchopulmonary segments varies from 18 to 20. This 239 causes the structure of central airways to vary among various people. To find the most common 240 structure, 50 CT images of various people were studied [37,38]. We found that the most common 241 structure is a lung with 18 segmental bronchi as shown in Fig. 1. Table 2 shows the name of each 242 segment fed by segmental bronchus indicated in Fig. 1.  Central airways shown in Fig. 1 can be generated using parametric L-system with seven rewriting 249 rules. Because the central structure is slightly different from the one reported by Davoodi and 250 Boozarjomehry [22], their production rules are modified to produce the central structure shown in 251 ( 3.2) where , , δ, β, and α are branch length, diameter, orientation angle and ratio of diameter and 253 length of daughter to its parent, respectively. +( ) and \( ) denotes orientation in azimuthal and 254 polar direction. If δ is positive (negative) the rotation will be counterclockwise (clockwise). All parameters are obtained by the analysis of CT images taken from Exact09 website [37]. 263 This dataset contains 40 CT scan images. The datasets used in this study are case01 and case37 264 which belong to adult male subjects. The airway segmentation module of 3DSlicer software was 265 used to extract central airways from CT images. The distance between two bifurcation points was 266 considered for branch length and average diameter of branch used as diameter. Polar and azimuthal 267 angles were estimated using the given branch orientation. All segmental bronchi (last branches of 268 central airways) are used as the axiom of the second section. 269

Dichotomous structure 270
To generate the dichotomous tree, parametric L-System is used. Because the rules proposed by 271 Davoodi and Boozarjomehry [22]are based on the bifurcating concept, we used their rules given 272 in Eq. (3.1). In each generation, a parent airway divides into two daughter airways in a single 273 branching plane. It should be noted that these rules are modified to include any rotation angle, 274 which is the angle between two successive branching planes. In addition, the process of finding 275 the parameters is different. θ and φ are branching angle, the angle between the branch and its parent, and rotation angle, 277 respectively. and are the normal vectors of the parents' plane and daughters' plane. The 278 schematic of two successive branching planes is shown in Fig. 2. and are the length and 279 diameter of the parent branch. is the ratio of branch diameter to that of its parent. 280 Turtle interpretation of the above rules is described in Table 3  Calculate the normal vector of daughters' plane and turn parent plane by angle φ −( ) Turn by angle in the clockwise direction +( ) Turn by angle in the counter-clockwise direction C Its value determines if the branching process continues from the daughter(s) branch or not 1 action of characters {A, E} is similar to F but the branching process terminates in A.

299
In the human bronchial tree, each branch divides into two branches with a reduction in diameter. 300 Mauroy et.al [32] showed that the best structure minimizing total viscous dissipation is a fractal 301 whose dimension is 3. In this ideal tree size ratio, the ratio of the daughters' dimensions to the 302 corresponding dimensions of their parent, is (1/2)1/3 ~0.79. In addition, they demonstrated that a 303 physically optimum tree does not have physiological robustness of the human bronchial tree. The 304 human bronchial tree has an average size ratio of 0.85, which is larger than the physically optimum 305 value, and with such a size ratio, the volume of the human bronchial tree is too large and its 306 resistance is too small. Some previous research aimed to design a bronchial tree based on 307 minimizing the volume and total resistance [41] or bifurcating resistance [16,22]. However, 308 physical optimization of the bronchial tree may be dangerous [32]. On the other hand, it is 309 necessary to model the bronchial tree to satisfy the physiological constraints. In what follows, we 310 suggest an intelligent procedure for finding the parameters of L-System to fulfill the physiological 311 and anatomical requirements of dimensions and orientations. 312 Parameters γ 1 and γ 2 must be less than 1 since the daughters' diameters are always smaller than 313 their parents. Murray[17] showed that Q=Cd3 minimizes the power dissipation in an ideal tree 314 assuming laminar flow. Uylings [42] extended the equation to Q=Cdn where 'n' is diameter 315 exponent and depends on flow regime, n=3 for laminar and n=2.34 for the turbulent regime. 316 Assuming a flow ratio 'r' the parameters γ 1 and γ 2 are: 317 (3.4) γ 1 = r 1 n γ 2 = (1 − r) 1 n Kitaoka et.al [16] found that n=2.8 leading to the best structure. Flow is not turbulent in all parts 318 of the dichotomous structure, so the value should be greater than 2.33. On the other hand, the 319 viscous forces are only dominant in the periphery of conducting airways thus the value should be 320 less than 3. 'r' is flow dividing ratio defined as the ratio of smaller flow in the daughters to that in 321 the parent. In this study, it varies from 0.15 to 0.5 indicating the structural asymmetry. As this 322 parameter gets closer to 0.5 the structure gets more symmetric. Since the human lung structure 323 becomes more symmetric with the generation, the parameter approaches to 0.5 with increasing 324 generation. Hence, this parameter is calculated by a normal distribution at each generation in such 325 a way that it remains between rmin and rmax, where rmin and rmax are functions of generation and 326 change after some generations to propel the structure toward symmetrizing. 327 To find the branch length an initial length (Lm > d) is considered. It is then examined whether 328 this branch intersects the Segment boundary and other branches during its growth or not. In 329 addition, this branch should be sufficiently far away from the Segment boundary so that there is 330 enough space for the growth of other conducting airways and the acinar region. Therefore, the 331 parameter D is defined as the minimum of the distance between the branch and other branches and 332 a fraction (rd) of its distance to the boundary. The latter was used to ensure the availability of 333 enough space for growth of acinar region which was ignored in the model of Davoodi and 334 Boozarjomehry[22] . Afterward, the branch length is adjusted in a value less than the minimum 335 of D and Lm according to Eq. (3.5). Kiatoka [16] suggested that the length to distance ratio ranges 336 from 0.17 to 0.34 by examining many computed tomography images. Because the length is 337 reduced using Eq. (3.5), the value of rd is set to 0.5. 338 where ω and is a measure of length reduction and is set at 0.04. Davoodi  the evaluation of branch distance to the boundary, the surface mesh became coarse as much as 359 possible. For this purpose, the triangular cells of different sizes have been used to keep the shape 360 of the segments unchanged as, for instance, indicated in Fig. 3. It should be noted that this 361 algorithm has been also used to obtain the distance between each branch with its surrounding 362 branches. For this purpose, each branch is surrounded by a cuboid which can be described by 12 363 triangular cells. Therefore, the distance can be approximated by the minimum distance of the 364 growing branch to the triangular cells of other branches. In addition, before the calculation of the 365 intersection of the branch with the other ones, the probability of its collision with them is examined. 366 In other words, for cases where the collision probability is not zero, the ray-triangle intersection 367 algorithm has been used. The intersection probability will be zero if the distance between the 368 centers of the branch and another branch is less than half of the sum of the length and diameter of 369 both branches. This makes the algorithm 20 times faster. It is clear that the airway length depends on the branch orientation that is determined after 374 finding the rotation angle of φ and the branching angle of θ. It should be noted that the parent and 375 its daughters are on the same plane. The schematic of parents and daughters have been shown in 376 Fig. 2. Accordingly, a normal vector to daughters' plane must satisfy the following equations. 377 where P, N, and N p are branch orientation of parent, normal vector to daughters' plane, and a 378 normal vector to the parents' plane, respectively. By solving the above equations, the elements of 379 N are obtained by: 380 Now, branch orientation can be determined by solving Eq. In this study, the rotation angle is fixed at 90°, which is consistent with the previous studies 386 [16,20,22]. The branching angle (θ) has a significant effect on branch length. 387 If an inappropriate value for θ is selected, many airways with small lengths are produced and 388 this prevents the growth of other airways and consequently, the development of conducting 389 airways would be stopped before reaching the acinar region. Besides, this causes the terminal 390 bronchioles with an inappropriate length to be generated. For example, an improper branching 391 pattern may lead to a bronchiole with a large diameter lies near the boundary, while its diameter 392 is not in the reasonable range. Furthermore, since it is in the vicinity of the boundary it cannot be 393 long enough and this causes too low length-to-diameter ratio. Such a scenario results in 394 inappropriate short branches. In addition, this may similarly cause small diameters. Furthermore, 395 the orientation of the branch affects the tree weight, fluid friction, and energy dissipation. 396 Therefore, it is essential to determine this parameter properly. Kitaoka[16]  However, all these methods were solely based on ongoing bifurcation, and this does not grantee 400 the minimization of the objective function for the whole structure of the human lung. Although 401 these mathematical relations decrease the resistance in the current bifurcation, it may result in lots 402 of bifurcations leading to the severe increase of total resistance. On the other hand, the human lung 403 is not physically optimized and such models are physiologically dangerous [32]. Therefore, an 404 intelligent procedure is required to lead the branches to the corresponding boundary so that the 405 anatomical constraints of the human lung are satisfied. Thus, instead of using any of the above 406 relations(e.g., Eq. (3.11)), the growth of each branch is examined in the several angles to obtain 407 enough information about the surrounding then the best angle is chosen based on some heuristic 408 rules used to satisfy anatomical and physiological constraints. For this purpose, the angle changes 409 by increment ∆θ between and and an angle is randomly selected at each subinterval 410 ([ , + ∆ )) by a normal distribution with a mean of and a standard deviation of ∆ /2. After 411 that, branch orientation is calculated for each θ using Eq. (3.9) or (3.10). To find the best branch 412 orientation among these potential directions, the following rules are used. It should be noted that 413 each direction has its diameter, length, and distance to the corresponding boundary whose 414 calculation algorithm mentioned previously. 415

Rule 1: if ⁄ ≤ . then the probability of being TB in the next few generations (i.e. up to 416 three generations) is high, thus a direction with the minimum length is chosen. This rule 417
prevents the production of TB's with too small diameter because if the minimum length was not 418 chosen, then there would be numerous growth cycles (generations) from this branch and this would 419 lead to branches whose diameters are significantly lower than their corresponding anatomically 420 meaningful values. 421

Rule 2: if . < ⁄ ≤ . then the probability of being TB in the next few generations is 422 medium, hence a direction with the minimum distance to the corresponding boundary is chosen. 423
This rule leads branches to the corresponding bronchopulmonary segment boundary and prevents 424 growing procedure stopped in the middle of the segment. 425

Rule 3: if the length of the selected direction is too small ( ⁄ < ), this direction is replaced by 426
another one that satisfies this condition and has the least length. This rule prevents the 427 production of branches with too small aspect ratio. 428 where ρ is a measure of emptiness around a growing branch which is defined in Eq. (3.13). 434

Rule 4: if the probability of being TB in a few next generations is low (i.e., neither of the above 429 rules is fired), the best direction is the one with the highest score defined in Eq. (3.12). This rule
where N is the number of branches with non-zero collision probability around the growing branch. 435 is the distance between the center of the i'th branch and that of the j'th branch. The second term 436 will be dominant in the last generations because of max((gen − 15) , 1). 437 As mentioned above, these rules generally conduct branches to bronchopulmonary segments 438 boundaries and less crowded regions. These rules primarily may cause few branches to be 439 generated because they may direct the branches forward rather than distribute them. To solve this 440 problem, the initial range of the branching angle ( and ) should be determined based on 441 its relative distance to the corresponding boundary in the axial and lateral directions to provide the 442 branch with the necessary information about its vicinity. These parameters tell the branch that it 443 is in a round space or a long thin region. In addition, to ensure that the lung surface and the 444 segments' interfaces are fed at least by one branch, the adjacent bronchopulmonary segments are 445 also considered in the determination of the angle range. Hence, conducting airways in all segments 446 should be generated simultaneously. If a portion of the interface is not fed by branches in the 447 neighbor segment, the angle is set in such a way that the chance of getting it fed through the 448 growing branch increases. Accordingly, the ratio of the distance between a branch and its segment 449 interface to the segment diameter along the direction of branch orientation ( ), and the same 450 for the lateral direction ( ), and the congestion difference on both sides of the interface (∆β) 451 are applied to obtain the range of the branching angle. ∆β is defined in Eq. (3.14) and calculated 452 for each cell of the interface. 453 where Acell and Ai are cross-section area of a triangular boundary cell and its adjacent airways, 454 respectively. 455 Although we understand that the range of branching angle depends on the mentioned parameters, 456 the exact relation is unknown. However, we know qualitatively the effect of each parameter on the 457 branching angle. To find the branching angle based on non-numerical information, the fuzzy 458 inference system is employed. The fuzzy inference system is a type of inference systems in which 459 the reasoning is done based on the facts and concepts whose truths are subjective (and /or cannot 460 be validated firmly). The logic behind such reasoning differs from conventional logic. In Fuzzy 461 logic, a fact might be partly true or false, while this is also the case for a rule or logical expressions. 462 The truth values of rules, logical expressions, and facts vary from 0 to 1. The design of a fuzzy 463 inference systems includes three steps: 1) Specifying the Fuzzy sets comprising the universe of 464 discourse of each variable (either input or output) used in The cascade fuzzy inference system is used to find the range of branching angles based on 472 previously mentioned affecting parameters. Figure 4 shows a schematic of the fuzzy system. 473  In the human bronchial tree, there are some branches called "links" whose orientations are in 495 the direction of their parent orientation. In order to include the links in the structure, the best 496 direction of each daughter is compared with the link and replaced via the following procedure. It 497 should be noted that only one daughter can be replaced and the procedure determines which one 498 is compared with the link. It worth noting that the "links" have been ignored in model of Davoodi 499 and Bozorgmehry[22]. 500 i) If 1 ( 2 ) =̂ and 2 ( 1 ) >̂, the daughter 1(2) is compared with the link to prevent 501 a very small angle between daughters. 502 ii) If θ 1 & θ 2 > θ min , the daughter with smaller / is compared with the link based on previous 503 rules. To ensure that the branches are well distributed, the generation of more than three subsequent 504 links are avoided. 505 According to the production rules, each branch is divided into two daughter branches provided 506 that the parameters of the daughters are determined. The orientation and dimensions of the 507 daughters are found via the following steps: 508 1. A random value is chosen for 'r' from a random number generator with a uniform 509 distribution. 510 2. The diameter of each branch is calculated using the diameter ratio given in Eq. (3.4). 511 3. Using a predefined value for rotation angle, φ, the normal vector of daughters' plan 512 calculated by Eq. (3.7). 513 4. The lateral and axial distances of each branch to the host boundary are calculated. 514 5. The value of parameter ∆β is obtained using Eq. (3.14). 515 6. The range of branching angle (θ min and θ max ) are determined through the proposed fuzzy 516 scheme. 517 7. Several branching angles are selected in the range of (θ min and θ max ), and the branch 518 orientation corresponding to each branching angle is calculated using Eq. (3.9) and Eq. 519 (3.10). 520 8. The distance to the boundary and parameter D is calculated for each potential branch 521 orientation. Then the corresponding branch length is calculated by Eq. (3.5). 522 9. The best orientation of each branch is selected among the potential orientations through the 523 proposed predictive rules. 524 10. At most one daughter branch is compared with the link and the best direction is chosen 525 using the same predictive rules. 526 11. The daughters are compared and if the branch with a larger diameter does not have a greater 527 length, the procedure is repeated from step 1. Because the branch with greater length has a 528 greater distance to the boundary and is located in the less crowded region, it should be 529 thicker to deliver air to the larger volume. 530 12. Although the predictive rules are designed to find the best direction, there may be few 531 terminal bronchioles with too small length. To remove such branches in the structure, the 532 algorithm goes back and reduces the diameter ratio in order for the paths to reach to terminal 533 bronchioles earlier and stop their growths in this direction before too short branches are 534 produced. 535 All parameters of the proposed procedure except

Results 543
In this section, the segmentation results of bronchopulmonary segments are presented. Then the 544 generated bronchial tree using the proposed method is evaluated through comparison of 545 geometrical properties with morphological measurements and the other alternate methods. The 546 analyses are performed for both central airways and the whole structure. In what follows, the 547 daughter branches with the smaller and larger diameters are called the minor and major branches, 548 respectively. The subscripts 'min', 'maj', and 'p' denote minor, major, and parent branches, 549 respectively. The results correspond to a total lung capacity of 5500ml. 550 The central airways have been generated via three generations of applying production rules. 551 The mean diameter of segmental bronchi, lobar bronchi, and total airways are 4.32mm, 7.8mm, 552 and 6.59mm which is comparable with the value of 4.61, 6.73, and 6.29 measured by Horsfield 553 et.al [13]. The insignificant difference between the values may originate from the difference in the 554 method of measurement and size of the lung. The geometric properties of central airways are 555 summarized in Table 6 Horsfield order are shown in Fig. 9. This figure indicates that the length doesn't decrease as 564 uniformly as diameter. The trend reveals that the length reduction has almost a linear trend for 565 small Horsfield orders. 566  The generation number is a measure of the distance of a branch and the trachea. The greater the 593 generation number reflects that the branches are further from trachea. Table 7 indicates the mean, 594 minimum and maximum generation number of terminal bronchioles in the lungs and each lobe. 595 The acini are appeared in the mean generation of 16.7, with a minimum of 8 and maximum of 25. 596 The number of TB's in the structure obtained by the proposed method is 30201 and is comparable 597 with those reported in the literature as shown in Table 7. The number of TB's obtained by  lies in the range of 26000 to 32000. Generation number of TB's is a 599 measure of the number of bifurcations from trachea at which an acinus may appear. Lobes that are 600 further away from trachea have a larger mean generation number of TB's as shown in Table 7. For 601 instance, the right lower lobe has the largest generation number because it has the furthest distance 602 from the trachea and is the largest lobe. 603 Although using generation number is useful for locating the branches in relation to the trachea, 606 it is not appropriate for asymmetric trees because this classification locates very dissimilar 607 branches in one generation. Therefore, the branches should be classified by order instead of 608 generation. Strahler and Horsfield orders are the most common classifying methods for 609 asymmetric trees. The number of branches is decreased with the increase in Strahler and Horsfield 610 orders and hence trachea has the highest Strahler and Horsfield orders. The relative increase in 611 each order is called branching ratio, Rb, and is defined as the antilog of the slope of the line 612 representing the number of branches versus the branching order in a logarithmic plot. The diameter 613 ratio (Rd) and length ratio (Rl) are calculated similarly and are measures of the reductions in 614 diameter and length, respectively. Figures 16-17 show the logarithmic plot of the number of 615 branches, diameter, and length against Horsfield and Strahler orders. Table 8 compares Rb,Rd,616 and Rl obtained by the proposed method and their corresponding published values. The value in 617 parenthesis is R squared (coefficient of determination) of the corresponding variable. Prefixes "S" 618 and "H" indicate the parameters calculated from Strahler and Horsfield ordering, respectively. 619 HRb and SRb are equal to 2 for symmetric trees. SRb becomes greater than 2 and HRb approaches 620 to 1 with an increase in asymmetry. SRb of 2.76 obtained in this study is higher than those given 621 by [18,21,47] and less than those predicted by [19,22]. HRb of 1.52 is consistent with the 622 published values. The difference in the reported values may be originated from the degree of 623 asymmetry of the methods and the shape of the lungs. Round lungs have a lower Rb than long thin 624 ones. Figure 15 shows a more linear trend comparing to those that exist in Fig. 16 and hence SRb 625 is more consistent throughout the lungs. However, every branch is not separately ordered in 626 Strahler ordering method. This makes a simpler description with loss of details compared to the 627 Horsfield ordering method in which each branch is separately ordered. To minimize resistance or 628 entropy production, the value of Rd is expected to be almost equal to Rb1/3[48]. The predicted 629 values for Rb and Rd well follow this relation (1.12% relative deviation for Strahler ordering and 630 0.25% for Horsfield ordering). 631  Because the diameters of three branches involving in a bifurcation are different, it is of interest 640 to observe the relation between them. In this study, the diameter decreased based on Eq. (3.4) 641 reflecting the relationship between flow rate and branch diameter in a living organ to make it 642 capable of overcoming the energy lost due to fluid friction in an optimal manner. For typical human 643 bronchial tree the diameter exponent, 'n', is approximately 2.92 [48], which is closer to that in the 644 laminar flow. We examined several values for 'n' in the range of 2.4 to 3.5 and found that the best 645 structure could be obtained with a value of 2.8 as proposed by Kitaoka et al[16] and used by 646 Davoodi and Bozorgmehry[22]. Weibel [29] showed that D/Dp is higher for generations >10 with 647 a minimum average value of 0.84. Therefore, the minimum value of 0.7 was considered for D/Dp 648 in generations>10. The stochastic behavior of the proposed procedure originates from the random 649 flow-dividing ratio (r) building slightly different structures with the same structural properties. It 650 is randomly set for each branch in such a way that the larger host volume is fed by the major 651 daughter r. D/Dp, Dmin/Dp, Dmaj/Dp, Dmin/Dmaj are consistent with those reported in the previous 652 studies. The obtained value of D/Dp is 0.856 which is higher than that predicted by the previous 653 methods [19][20][21][22] and slightly less than that measured by Weibel [29]. In an ideal fractal tree, this 654 ratio is 0.79 (21/3). However, Mauroy et al [32]showed that the human bronchial tree is not 655 physically optimal and the ratio in the real tree is about 0.85 which is approximately equal to that 656 obtained in this study. Diameter reduction in the human bronchial tree is greater than that predicted 657 by the previous models and that is due to the fact that the physiological robustness of human lung 658 is obtained at the expense of its deviation from the physically optimum structure. This 659 physiological robustness makes the lungs functional in the case of slight damage in its structure. 660 That is, the diameter ratio in the proposed method is greater than that in the previous models. 661 Furthermore, D/Dp predicted by the proposed method has been compared with those measured by 662 Weibel[29] and Yeh and Schum[11] in Fig. 18. This figure confirms the ability of the proposed 663 procedure in the prediction of real data.
/ calculated by the present model is much smaller 664 than those reported by Bordas et.al [20]  The predicted mean value of the length-to-diameter ratio is inline with those calculated in the 672 previous studies as indicated in Table 8. The higher standard deviation (SD) may be originated 673 from considering bronchopulmonary segments as host volume while the tree was grown in the 674 lungs or lobes in the previous studies. For instance, the inferior lingular is a long thin segment. 675 Therefore, airways along axial direction should have a large length, whereas those in lateral 676 direction have smaller length, this results in greater standard deviation compared to the previous 677 models. However, the coefficient of variation is comparable to what estimated by Tawhai et.al 678 [19]. 679 / has a higher standard deviation than / because the length doesn't decrease as uniformly 680 as diameter. Figures 16-17 Davoodi and 687 Bozorgemhry [22] did not report the value of DTB, we obtained too small TB's using their method. 688 This is a major shortcoming because the small dimensions of TB's consequently affect the size of 689 the acini fed by them and subsequently underestimate the lung volume and gas exchange area too 690 much. 691 The TB length (LTB) of 0.87±0.40mm has been predicted by the proposed method. The LTB of 692 0.819±0.14mm, 0.8±0.354mm, and 0.93±0.43mm reported by [39], [23], and [40] using the 693 direct measurements on TB's. However, the values of 1.34±0.44mm and ~1.8mm predicted by 694 [20] and most of the previous models [19,31,50]. The LTB predicted by the present method is less 695 than those reported by the previous models and is more consistent with the experimental data. This 696 is because, in the present method, the airways grow inside the bronchopulmonary segments which 697 restrict the growth and length of TB's more than what lobes or lungs do. 698 The mean branching angle in our model is 41.4°. The mean branching angle is about 47° based 699 on the data of 20 lung casts [31] and 41.26° reported by [41] which is close to that calculated by 700 the proposed method. The calculated SD based on the data of [41] is 23.26°, which is almost the 701 same as what obtained in this study. Due to the asymmetric behavior of the tracheobronchial tree, 702 all proposed models experience relatively large standard deviations except the model of [20] with 703 almost no SD. Although the model of [20] has an asymmetric manner, it seems that the asymmetric 704 pattern is almost identical in all bifurcations as obvious from reported SD for other properties. 705 The calculated volume of the tracheobronchial tree is 133ml. The volume of 145ml was measured 706 by Weibel et.al [29]. Kitaoka et.al [16]and Ismail et.al [51]estimated the volume of 175ml and 707 132.9ml, respectively. The difference between the reported values may correspond to the fact that 708 they are obtained based on the lungs of different subjects. The anatomical dead space measured by 709 the Fowler method is 156±28ml. excluding the mouth volume of 40ml (estimated as 70% of the 710 maximum volume of mouth[52]), the volume of conducting airways is expected to be 116±28ml. 711 Therefore, the value obtained by the proposed method is reasonable. The percentage volume of 712 airways in each lobe is shown in Table 9. This table indicates that the volume of conducting 713 airways is related to the volume of the corresponding lobe. 714 715

Discussion 719
In this study, to generate the human tracheobronchial tree, we tried to provide a robust rule-based 720 method whose parameters were determined using an intelligent procedure. The bronchial tree was 721 produced such that its various parts fit into bronchopulmonary segments obtained by analysis of 722 CT images and anatomical slices. Various parts of the tree obtained by the proposed algorithm are 723 completely contained by their corresponding host volumes in a homogenous manner while 724 fulfilling various statistical indices reported by anatomical and physiological research articles. 725 The spatial hindrance during growth of an airway due to the presence of surrounding airways 726 have been included in three-dimensional manner in this study, whereas it corresponding 727 counterpart was handled in one-dimensional in the previous works. Furthermore, the proposed 728 method predicted the more accurate dimensions of TBs and D/Dp ratio than their counterparts 729 obtained by previous asymmetric models. 730 In the current study, the TB's are distributed near the bronchopulmonary segments, which is 731 what happens in reality. Longitudinal pathways, the distances along the airways from the trachea 732 to TB's, have a great contribution to the resistance against the fluid flow and subsequently on the 733 ventilation and concentration distribution. Because TB's have been located near the boundaries of 734 lungs or lobes in the previous studies, the values of Longitudinal pathways should have a 735 significant deviation from reality but the values have not been reported. 736