Developmental model of an automatic production of the human bronchial tree based on L-system

computer methods and programs in biomedicine 132 (2016) 1–10

j o u r n a l h o m e p a g e : w w w. i n t l . e l s e v i e r h e a l t h . c o m / j o u r n a l s / c m p b

Developmental model of an automatic production of the human bronchial tree based on L-system Amirabbas Davoodi, Ramin Bozorgmehry Boozarjomehry * Chemical and Petroleum Engineering Department, Sharif University of Technology, Azadi Av., Tehran, Iran

A R T I C L E

I N F O

A B S T R A C T

Article history:

Background and objective: The human lungs exchange air with the external environment via

Received 10 July 2015

the conducting airways. The application of an anatomically accurate model of the conduct-

Received in revised form

ing airways can be helpful for simulating gas exchange and fluid distribution throughout

16 March 2016

the bronchial tree in the lung.

Accepted 19 April 2016

Methods: In the current study, Lindenmayer system (L-system) has been formulated to generate the bronchial tree structure in a human lung. It has been considered that the structure

Keywords:

of the bronchial tree is divided into two main segments: 1) The central airways (from the

Bronchial tree

trachea to segmental bronchi) and 2) the dichotomous structure (from segmental bronchi

L-system

to terminal bronchioles). Two sets of parametric rewriting rules which can be used to develop

Lung

central and peripheral airways have been proposed; the first set used to develop central airways

Rewriting rule

consists of seven rules, while the second rule set contains four rules. Results: The proposed model is capable of generating bronchial tree inside the volume of the host lung; and comparison of the resulting model with those reported in the literature shows that the morphometric characteristics of L-system structure are in good agreement with their corresponding experimental data. Conclusion: The resulting model can be used to obtain a mathematical model required for the study of transport phenomena occurring in the lung during respiration. © 2016 Elsevier Ireland Ltd. All rights reserved.

1.

Introduction

The most important role of the lungs is exchanging air between organisms and environment. The lungs transport air via the bronchial tree (from the trachea to the terminal bronchioles) and the gas exchange with the blood takes place in the alveoli air sacs [1,2]. The bronchial tree in the human lungs and those of other mammals is a complex branching structure. The human bronchial tree is an asymmetric branching structure with the tendency to form dichotomous division into two daughter

airways with different dimension, usually diameters and branching angles [3]. Structural models of the bronchial tree provide useful results for investigating the correlation between structure and pulmonary function. Weibel’s symmetric structural airway model is one of the most well-known airway models that contain airway dimensions and connectivity [4]. This model provides a good basis for studying lung functions and gas exchange; however, it obviously cannot be appropriate for investigation of the effect of anatomical asymmetry and spatial location of airways on transport problems. Horsfield and Cumming have studied the

* Corresponding author. Chemical and Petroleum Engineering Department, Sharif University of Technology, Azadi Av., Tehran, Iran. Tel.: +98 21 6616 6445; fax: +98 21 66022853. E-mail address: [email protected] (R.B. Boozarjomehry). http://dx.doi.org/10.1016/j.cmpb.2016.04.021 0169-2607/© 2016 Elsevier Ireland Ltd. All rights reserved.

2

computer methods and programs in biomedicine 132 (2016) 1–10

functional effect of asymmetry of the bronchial tree by preparing a cast with thermosetting resin [5]. They have demonstrated that numbering orders upward starting with the end branches is more suitable than counting generations downward from trachea for asymmetrical structures. They proposed two ordering classifications, namely Horsefield and Strahler. According to the Horsfield ordering method, the terminal branches are numbered the order one and at each branch junction, the parent branch gets one order more than the highest order of the related daughter branches [5]. Based on the Strahler ordering [5], the parent branch is assigned one order higher than two daughter branches of the same order; otherwise the parent branch continues with the same order as the daughter branch with the highest order. Next, contiguous branches of the same order are then replaced by a single branch. Horsfield et al. [6] have described a conducting airway model of the bronchial tree with regular asymmetry based on delta, the difference between the order of the daughter branches. Horsfield and Thurlbeck [7] developed a computer program to simulate the human bronchial tree. They ordered trees by the methods of Horsfield and Strahler and used the aforementioned delta model to investigate the effect of some morphological parameters on the form of trees. There are some similar attempts to generate geometric airway models, but these models are still limited to two-dimensional structures [8–10]. Kitaoka et al. [11] proposed a deterministic algorithm which used nine basic rules and four complementary rules to generate a three-dimensional model of the human airway tree that can supply fluid within the organ. Tawhai et al. [12] presented a method for generating three-dimensional bronchial airway which grows into finite element volume mesh of the five lobes of the human lungs. Tawhai et al. [13] derived the models of the human and sheep bronchial trees directly from the detailed anatomic data obtained from multidetector row X-ray-CT (MDCT) scans and used Tawhaiet al. [12] model to generate anatomically based airway models that start from the periphery of the MDCT-segmented airways and fill the volume meshes. Because of the rapid advancement in multislice spiral CT scanner technology, several methods have been developed for segmenting and analyzing the 3D airway trees [14–18]. The evaluation of these methods shows that the precision of the CT images plays an important role on the obtained structures of the bronchial tree; besides, the results do not include any concept corresponding to the growth of the airway branches. Three-dimensional models of the bronchial tree have been widely used for various purposes, from CFD analysis of airflow and aerosol deposition in human lung [19–21] to prediction and assessment of airway functions [22] and tracing airway paths for virtual bronchoscopy [13,23–26]. All in all, the application of an appropriate model with a high level of anatomic detail is essential for the study of the aforementioned fields. Lindenmayer system or L-system, introduced by Aristid Lindenmayer [27,28] is a simulation approach for the development of the multicellular organisms. L-systems can be classified as rule-based technique that uses a set of grammar or generative rules to create branching structures. Originally, L-system was proposed as a mathematical theory of plant development, and it did not cover geometric aspects [27,28].

Subsequently, several geometric interpretations of L-systems were proposed with a view to turning them into a mature framework that could model plant shape and growth [29]. To the best of our knowledge, although L-systems have been applied to create the branching structure of various systems like botanical plants [29–34] and blood arteries [35,36], no similar attempts have been carried out using this method in order to generate distributed models of the bronchial tree. This technique makes it possible to define some characteristics of a branch, e.g. diameter, length, and branching angle, which can be calculated in parallel and with independent equations. This approach provides a framework by which the airways in lung with structural abnormalities can be constructed. Furthermore, it can address the evolution of the constructed airway structure based on the growth of the lung volume which happens throughout the time (from infantry to maturity). In the present work, a novel method in order to generate bronchial trees growing into any given lung space is developed. The appropriate L-systems with some rewriting rules are introduced to generate human conducting airway. The proposed algorithm is investigated for the assessment of the applicability of this method in order to generate the human bronchial tree.

2.

Theory

2.1.

Basic L-systems

L-systems are parallel rewriting systems, in which complex strings are defined from initial string, called the axiom by using a set of production rules. Each rule includes a predecessor and a successor by which the predecessor is iteratively replaced by the successor in the resulting string. If no production rule is explicitly specified for a given predecessor, the latter is reiterated in the newly generated string. Assigning geometrical interpretations to different symbols, for instance based on turtle interpretation, lead to the interpretation of the final string as a spatial object [29]. The simplest class of L-systems is called DOL-systems (deterministic and context-free L-system). Formal OL-system (context-free L-system) is an ordered triplet G = V , ω , P where V is an alphabet, if the set of all words and all nonempty words over V is denoted by V* and V + , respectively, then ω ∈ V + is the axiom and P ∈ V × V * is a finite set of production rules. The production rule is written as the following format:

P : predecessor → successor

2.2.

(1)

Parametric OL-system

A parametric OL-system is an ordered quadruplet G = V , Σ, ω , P where V is the alphabet of the system, Σ is a set of formal parameters, ω is the axiom, and P is a finite set of production rules. The parametric production rule has the following format:

P : predecessor ( parameters) : condition → successor ( parameters) (2) The symbol ‘:’ is used to separate the predecessor and the condition.

computer methods and programs in biomedicine 132 (2016) 1–10

2.3.

3

Bracketed OL-system

Lindenmayer used strings with bracket to represent the axial tree [27,28,37]. An axial tree from alphabet V can be represented by string of symbols over alphabet VB = V ∪ {[,]} . The brackets are used to restrict the branches. Further details on all types of L-system with several examples can be found in [29].

3.

Method

3.1. Application of L-system in conducting airway modeling The human bronchial tree is a geometric structure, which is propagated within two lungs to supply air for respiratory units. Moreover, the generated model should consider the anatomical constraints of structure dimensions and orientation. Therefore, the bronchial tree has been divided into two main segments, and the parametric OL-system, constructed by the physiological laws of the airway branching, is used to generate each segment.1 The branch dimensions are defined by means of the length and the diameter, while the branch orientations are calculated through both polar and azimuthal angles. The first section, called central airways, starts from the trachea and includes main bronchus, lobar bronchi, and segmental bronchi. The second section consists of dichotomous branches which starts from the periphery of central airways and stops at the terminal bronchioles. All branches with similar lobar bronchi grow into their corresponding lobe. In each lobe, the process of generating dichotomous structure from segmental bronchi occurs simultaneously; however, each lobe demands its own separate computational effort.

3.2.

Derivation of L-system for the central airways

The central airways are generated by applying the parametric OL-system method, including seven rewriting rules. The schematic structure of the central airways, ordered by the Strahler method has been shown in Fig. 1. This structure can be formulated by the following production rules which produce branches with lower order in each generation:

ω : T (L0, D0 ) p1 : T ( L, D) → + (δ ) \ (δ ) E (α L, β D) [ + (δ ) \ (δ ) B (α L, β D) + (δ ) \ (δ ) C (α L, β D)][ + (δ ) \ (δ ) G (α L, β D) + (δ ) \ (δ ) M (α L, β D)] p2 : B ( L, D) → E (α L, β D) [ + (δ ) \ (δ ) K (α L, β D) + (δ ) \ (δ ) M (α L, β D)] p3 : C ( L, D) → E (α L, β D) [ + (δ ) \ (δ ) M (α L, β D)][ + (δ ) \ (δ ) K (α L, β D) [ + (δ ) \ (δ ) K (α L, β D) + (δ ) \ (δ ) N (α L, β D)] p4 : G ( L, D) → E (α L, β D) [ + (δ ) \ (δ ) K (α L, β D) + (δ ) \ (δ ) N (α L, β D)] p5 : M ( L, D) → E (α L, β D) [ + (δ ) \ (δ ) N (α L, β D)][ + (δ ) \ (δ ) M (α L, β D)] p6 : K ( L, D) → E (α L, β D) [ + (δ ) \ (δ ) A (α L, β D)] p7 : N ( L, D) → E (α L, β D) [ + (δ ) \ (δ ) A (α L, β D)][ + (δ ) \ (δ ) (3) A (α L, β D)][ + (δ ) \ (δ ) A (α L, β D)] 1 The code and software tutorial can be found in the Supplementary Data.

Fig. 1 – The schematic structure of central airway ordered with Strahler method.

where L0 and D0 are the length (mm) and diameter (mm) of the trachea, respectively, δ denotes the orientation angle, α and β are the ratio of daughter to parent length and diameter, respectively. The trachea (T(L0,D0)) has been considered to be the axiom, where the final string is constructed through three generation of applying central airways production rules. All parameters can be measured manually from chest Computed Tomography (CT) images. The images were taken from EXACT‘09 website [38], the publicly available dataset of 40 scans, including 20 scans for training and 20 scans for testing. The images named CASE01 through CASE40. A detail description of these cases was presented in [18]. The only dataset used in this study corresponds to CASE01 of these datasets. The airway segmentation module [39] was used to segment central airway from CT images by Slicer4 [40]. Then, the segmented tree was applied to obtain anatomic characteristic measurements in which the branch point regions were neglected, and the cylindrical parts of the corresponding segmented airway were considered. The distance between each bifurcation point was estimated for the respective branch length, and the average value of minor and major diameter was used as the branch diameter. In addition, the angle between the polar and azimuthal direction of the parent and daughter branches was estimated. The obtained data are sequentially assigned to the parameters of rewriting rules in each generation. The geometries of all segmental bronchi are stored as properties of the axiom for the next section.

3.3.

Designing dichotomous structure

In the following, four parametric OL-systems based on bifurcating concept are proposed to produce dichotomous tree, whereby in each generation, a segment airway branch divides into two branches in a single plane then each branch continues the same, and so on. It is assumed that the rotation angle, the angle between the two successive branching planes, is 90°

4

computer methods and programs in biomedicine 132 (2016) 1–10

Table 1 – States of parameter π. State 0 1

Action The angle varies in the polar direction The angle varies in the azimuthal direction

which is consistent with the value proposed in previous studies [13].

ω : F (l0, d0, 0, 0) p1 : F (l, d, π , 0) → E (l, d) [ − (θ1 ) F ( Γ ( D) , γ 1d, 1 − π , c)][ + (θ 2 ) F ( Γ ( D) , γ 2d, 1 − π , c)] p2 : F (l, d, π , 1) → E (l, d) [ − (θ1 ) A ( Γ ( D) , γ 1d)][ + (θ 2 ) F (Γ ( D) , (4) γ 2d, 1 − π , c)] p3 : F (l, d, π , 2) → E (l, d) [ − (θ1 ) F ( Γ ( D) , γ 1d, 1 − π , c)][ + (θ 2 ) A ( Γ ( D) , γ 2d)] p4 : F (l, d, π , 3) → E (l, d) [ − (θ1 ) A ( Γ ( D) , γ 1d)][ + (θ 2 ) A ( Γ ( D) , γ 2d)]

Fig. 2 – Branch bifurcating that happens in a single plane.

where, l and d denote the length (mm) and diameter (mm),respectively. θ is the branching angle, and γ is the ratio of branch diameter to the parent diameter. In addition, π and c are two logical parameters that determine the branching plane and branching process scenario, which are described in Tables 1 and 2, respectively. Furthermore, Γ (mm) is the function calculating the branch length based on D(mm), which is the minimum distance from starting point of the daughter branch to the lobe boundary or other branches along the growing direction.

where Q is the fluid flow rate, C is a constant depending on the organ and fluid, and n is the diameter exponent [43]. Uylings [44] argued that n depends on the flow regime, since in the central airways the flow is not laminar (it is either turbulent or transitional), n must be close to 2.3, whereas for laminar flow that exists in the periphery of the conducting airways, n must be close to 3.0. Kitaoka et al. [45] estimated n to be 2.8 by examining morphometric data of four mammalian airway trees published by Raabe et al. [46]. Hence, in this study n is assumed to be 2.8. The fluid flow is conserved before and after branching (Fig. 2), and the flow-dividing ratio (r) is defined as the ratio of the flow in daughter branch receiving the smaller flow to that in the parent, which ranges between 0 and 0.5 [11]. If r is close to the upper limit, then the generated structure will be more symmetric. The branching angle (θ ) and (γ ) can be derived as follow (Q 1 ≤ Q 2 )

Γ ( D) = min ( Γ max, D) × exp ( − ξ D)

γ1 = r

(5)

1

n 1

where Γ max (mm) is the maximum step for length of a branch, and ξ (mm) is the fractional factor that controls the branch length proportional to the parameter D. The parameter D is calculated in two stages: first, the distance between starting point of the branch and the lobe boundary is computed. The mathematical models for boundary of all five lobes are required to calculate distances of each branch to the corresponding boundary. Therefore, two lungs were segmented by the method presented in [40,41], and lobe fissures were segmented manually by Slicer4 [40]. Then, a set of points was selected on each segmented surface, and thin-plate smoothing spline method was applied to represent the lobe boundaries mathematically. In the next stage, it is assessed whether any branch exists in the path of the daughter branch along growing direction or not. Five points with the same distance on centerline of each branch have been considered as nominees of the existence of that branch in the daughter branch path. If there is any point in the daughter branch path, the distance from that point will be calculated. Obviously, the points of the ancestor branches are excluded, and are not checked in stage 2. Finally, the minimum value of those obtained in stages 1 and 2 is considered for parameter D. The calculation of parameters d, γ , and θ are based on a power law relationship between airway diameter and flow method, which was first proposed by Murray [42]

Q = Cdn

(6)

γ 2 = (1 − r )

(7) n

θ1 = cos−1 ([1 + γ 14 − γ 22 ] 2γ 12 )

θ 2 = cos−1 ([1 + γ 24 − γ 12 ] 2γ 22 )

(8)

It should be noted that, in this study, based on the convention used in this article, subscript 1 always refers to the daughter branch which is produced first from the successor string in the dichotomous branching structure rewriting rules. Kitaoka et al. [11] set up a threshold value for the lower limit of the flow-dividing ratio as the value of 0.05. In this study, r has been considered as a random variable that changes between 0.05 and 0.5 to add stochastic characteristic to the given algorithm. According to the set of rules developed for the dichotomous structure, each branch is divided into two branches in each generation, and the dimension and orientation of each branch is calculated through the following steps: 1) A random value is selected for the flow-dividing ratio (r). 2) The diameter and angle of each branch are calculated by Eq. 7 and Eq. 8, respectively. 3) According to the orientation of the two daughter branches, the parameter D is calculated for each branch. If the branch with larger diameter has higher value for D, the length of daughter branches will be calculated by Eq. 5; otherwise, the procedure is repeated from step 1 until the branch with

5

computer methods and programs in biomedicine 132 (2016) 1–10

Table 2 – States of parameter c. State 0

Branching process can be continued in both daughter branches Branching procedure should terminate in daughter branch number 1a Branching procedure should terminate in daughter branch number 2 Branching procedure should terminate in both daughter branches

1 2 3 a

Action

Table 3 – Symbols of commands in turtle interpretation for central airways rewriting rules. Command T (L, D) a + (δ )

\ (δ )

The branch number is based on the convention described previously. a

larger diameter reaches higher value in D. Since the distance from the lung wall is large, the branches have larger volume to occupy and deliver air, and hence the branches tend to be thicker. The parameter c determines the daughter branch in which the branching process can be continued, various values of this parameter and their corresponding effect on the branching mechanism are shown in Table 2. It should be noted that limits on the length and diameter dictate the value of this parameter at various growth cycles. The length of human terminal bronchiole ranges from 1 to 1.72 mm. A value of 1.2 mm is imposed for the limit on the branch length. Furthermore, the ratio of diameters of the daughter branch and the trachea is considered as a threshold for the lower limit of the branch diameter. A value of 0.035 is assigned for this threshold, which is estimated from the results presented in [13]. If the length and diameter of the daughter branch are smaller than or equal to the imposed limits of length and diameter, the daughter branch will be considered as terminal bronchioles and will be shown by the character {A} in the dichotomous production rule. The process continues until all characters {F} are replaced by {A} in the generated string.

3.4.

Action Move forward with a step of length L and draw a cylinder with length L and diameter D. Turn in the azimuthal direction by angle δ . If δ is positive, the turtle is turned in counterclockwise direction and if δ is negative, the turtle is turned in clockwise direction. Turn in the polar direction by angle δ . If δ is positive, the turtle is turned in counterclockwise direction and if δ is negative, the turtle is turned in clockwise direction.

Actions of characters {B,C,E,G,K,M,N} are similar to the action of character T.

4.1.

The central airway structure

Fig. 3 shows three-dimensional model of the central airway tree constructed using the central airways rewriting rules as well as the airways segmented from the CT images of the human male. The structure is generated in three generation based on

Table 4 – Symbols of commands in turtle interpretation for dichotomous branching structure rewriting rules. Command F (l, d) a + (θ ) − (θ ) a

Action Move forward with a step of length l and draw a cylinder with length l and diameter d Turn in counterclockwise direction by angle θ Turn in clockwise direction by angle θ

Actions of characters {A,E} are similar to that of character F.

Turtle interpretation

The generated L-system can now be used to model the bronchial tree. Turtle interpretation should be introduced to interpret the resulting string. The turtle interpretations of both central and dichotomous airway structures are implemented in spherical coordinate system. The command symbols used in production rules for central airways and dichotomous structure are listed in Tables 3 and 4, respectively.

4.

Results

The bronchial tree generated by L-system has been analyzed with morphometric studies and published mathematical models. Analyses have been done for the central and dichotomous structures separately. The daughter branches with larger and smaller diameters have been defined as major and minor branches, respectively.

Fig. 3 – Three-dimensional structure of the central airway reconstructed from X-ray-computed tomography cross-sectional image from the supine human subject (a) and generated by L-system (b).

6

computer methods and programs in biomedicine 132 (2016) 1–10

Table 5 – Geometric analysis of the human central airway structure.

Table 6 – Geometric analysis of the human dichotomous branching structure.

Parameters

Parameters The dichotomous branching model

θ (total) L/D D/Dparent Dminor/Dmajor

The central airway model

Published

37.38° ± 15.75° 3.42 ± 1.61 0.63 ± 0.16 0.69 ± 0.17

36.11° [13]; 39° [47] 3.04 [13]; 3.09 [47] 0.85 [13]; 0.82 [47] 0.71 [13]; 0.83 [47]

Results are mean ± SD. L/D, length-to-diameter ratio; θ, branching angle. Values for the same parameters from previous studies in the final column are followed by their corresponding references.

Strahler order. The mean diameters of total airways, lobar bronchi, and segmental bronchi are 4.8 mm, 6.3 mm, and 2.9 mm, respectively. The geometric parameters evaluated for the central airways are summarized in Table 5. The mean value of θ is approximately in the middle of the range of published values [13,47]. The mean θ calculated for branches with Strahler order 1, 2, and 3 is 35.8°, 39.07°, and 40.59°, respectively. The results show the increase of θ values with the increase of branch order. The mean value of length-to-diameter ratios (L/D) is slightly larger than what presented in previous studies [13,47]. The deviation of maximum value and minimum value of the estimated L/D from its corresponding mean value are smaller than what reported in [13]. The mean value of the ratios of D/Dparent through generations 1–3 is 0.7, 0.71, and 0.58, respectively. The mean D/Dparent decreases with the increase of generation number. However, the total mean value is smaller than its corresponding measured one from previous studies [13,47]. Classification of minor and major branches enforces the mean ratios of minor to the major diameter (Dminor/Dmajor) to be smaller than 1. The value obtained in this study is generally smaller than the range reported in previous studies [13,47].

4.2.

The dichotomous structure

Different values of Γ max and ξ have been used in order to generate dichotomous tree in each of the lobe volumes. According to the resulted data, the optimal values of Γ max and ξ have been found to be 5.5d (5.5 times the diameter) and 0.3, respectively. The geometric parameters obtained based on the parameters used to generate dichotomous structure are listed in Table 6. The mean value of θ is smaller than its corresponding published values [48]. The maximum and minimum values measured for θ are 69.14°, and 6.32°, respectively. This range results in large standard deviations for θ. The calculated mean value of length-to-diameter ratio (L/ D) for the dichotomous structure is larger than what was reported in previous studies [4,11]. Also, this deviation is observed for the mean value of length-to-diameter ratio of the minor and the major child branches. L/Dminor is closer to the total mean value and the values of L/D, L/Dminor, and L/Dmajor follow the trend of results (L/Dminor < L/D < L/Dmajor), which have been published in [13]. The mean values calculated for parameters D/Dparent, Dminor/ Dparent, Dmajor/Dparent, and Dminor/Dmajor are approximately consistent

θ (total) L/D L/Dminor child L/Dmajor child D/Dparent Dminor/Dparent Dmajor/Dparent Dminor/Dmajor L/Lparent

36.1° ± 14.88° 4.05 ± 0.41 4.04 ± 0.43 4.07 ± 0.39 0.76 ± 0.14 0.64 ± 0.15 0.88 ± 0.09 0.75 ± 0.19 0.78 ± 0.7

Published 37.28 [48] 2.8–3.25 [4]; 3.08 [11]; 2.92 [13] 2.88 [13] 2.96 [13] 0.79 [4,13]; 0.78 [47] 0.69 [13] 0.88 [13] 0.74 [47]; 0.81 [13] 0.81 [13]

Results are mean ± SD. L/D, length-to-diameter ratio; θ, branching angle. Values for the same parameters from previous studies in the final column are followed by their corresponding references.

with their corresponding values reported in previous studies [4,13,47]. The values calculated for Dminor/Dparent and Dmajor/Dparent demonstrate the flow-dividing ratio (r ) mainly varies between 0.28 and 0.31. The mean value of L/Lparent is smaller than the one which is consistent with the method used to calculate the branch length (Eq. 5). However, the value is smaller than what was measured in the previous study [13].

4.3.

The L-system bronchial tree model

The five dichotomous structures were merged with the central structure to produce a full bronchial tree model. The final model of the bronchial airway is shown in Fig. 4. Branches distal to lobar bronchi are shown with different colors: the right upper and left upper lobes are red, right middle lobe is green, and right lower and left lower lobes are blue. The average generation number from the trachea to the terminal bronchioles is 15.96 with a minimum value of 7 and maximum value of 23. Weibel [4] measured an average of 16

Fig. 4 – Anterior views of the bronchial airway model obtained by L-system.

computer methods and programs in biomedicine 132 (2016) 1–10

7

Fig. 5 – Log plots of the number of branches, mean length, and mean diameter against Strahler Order. The linear trendline of each plot is drawn as dashed line.

generations from the trachea to the terminal bronchiole, whereas Horsfield and Cumming [5] found an average of 14.6 generations from the trachea to the lobular branches. Kitaoka et al. [11] gave an average of 17.6 for their whole model. Tawhai et al. [12] reported the average number for each lobe whose mean is around 16. There are 29,973 terminal bronchioles in the generated model. Horsfield and Cumming [5] and Kitaoka et al. [11] estimated 27,992, and 27,306 terminal bronchioles, respectively, whereas Tawhai et al. [12] estimated 29,445. The result is in the range of 26,000–32,000, which have been estimated by HaefeliBleuer and Weibel [49]. Logarithmic plots of the number of branches, mean length, and mean diameter against Strahler order and Horsfield order are shown in Figs. 5 and 6, respectively. As shown in Figs. 5 and 6, the number of branches with the same order increases as the order decreases, whereas the mean length and mean diameter

increase as the order increases. The trachea has the highest order in the structure. It is at the order of 11 based on Strahler ordering scheme and at the order of 29 based on Horsfield ordering method. The antilog of the slope of log(number of branches) against the order is branching ratio, Rb. The length ratio Rl and the diameter ratio Rd are calculated similarly from the antilog of the slope of log(mean length) or log(mean diameter)against the order [50]. Tables 7 and 8 compare Rb, Rl, and Rd, which have been calculated from Strahler ordering and from Horsfield ordering of the model obtained by the propose method against their corresponding published values [5,12].

5.

Discussion

In the present study, a developmental method for the automatic production of the human bronchial tree based on the

Fig. 6 – Log plots of the number of branches, mean length, and mean diameter against Horsfield Order. The linear trendline of each plot is drawn as dashed line.

8

computer methods and programs in biomedicine 132 (2016) 1–10

Table 7 – Branching, length, and diameter ratios calculated using Strahler ordering. Ratio Rb Rd Rl Rb1/3

The L-system model

Published

2.81 1.51 1.47 1.41

2.8 [5]; 2.35 [12] 1.43 [5];1.32 [12] 1.4 [5]; 1.34 [12]

Values from previous studies in the final column are followed by their corresponding references.

Table 8 – Branching, length, and diameter ratios calculated using Horsfield ordering. Ratio Rb Rd Rl Rb1/3

The L-system model

Published

1.44 1.14 1.11 1.13

1.39 [12] 1.11 [12] 1.11 [12]

Values from previous studies in the final column are followed by their corresponding references.

L-systems has been proposed and implemented. The proposed method is a rule-based method that helps adding morphometric features to the airway tree structure during its development. The method contains a collection of rewriting rules whose predecessors in the string (obtained at the end of each growth cycle) are iteratively replaced by successors, and the structure is produced in the format of a complex string. Then, the string is interpreted, and a 3D image of the model is constructed through the turtle interpretation mechanism. As shown earlier, the human bronchial airway was divided into two segments; the central airway and the dichotomous structure. First, the central airway was constructed by 7 rewriting rules and then the dichotomous structure was generated by 4 rewriting rules in periphery of the central airway. The dichotomous structure consists of branches which have grown inside the lobes volumes. The rewriting rules of the proposed L-system were applied for each lobe separately. In contrast to the Weibel’s symmetric model and the Horsfield’s model [4,6], the proposed model is filling the volume of the host lung similar to what presented in [11,12]. The characteristic anatomic data of the structure of the complete bronchial tree, obtained by the proposed method, are in good agreement with their corresponding values reported in the literature. The dimensional parameters of the central airway were extracted from the thoracic human CT images for an individual. The geometry analysis shows (Table 5) only the mean value of branching angle is close to published values [13,47]. The deviation in other parameters (L/D, D/Dparent, Dminor/Dmajor) may originate from the differences exist in structures which have been evaluated. For example, the central airway constructed by the proposed method has 32 branches while the MDCT airway has 121 branches [13]. The generation of the dichotomous structure is based on the bifurcation; that is, there are no divisions into three (trichotomy) or more branches (polychotomy). Although the occasional trichotomy may be found in the lung due to

shortening of a branch to zero length [4,5], because of nonavailability of anatomic data for validation, it has been neglected in this study. Values for branching, length, and diameter ratios were calculated using both Strahler and Horsfield ordering. The evaluation shows that for minimal resistance or for minimal entropy production it is expected to have Rb1 3 = Rd = Rl [50]. In Strahler ordering, Rb is close to the Horsfield and Cumming [5] model, in reverse, it is larger than Tawhai et al. [12] results. The ratios for branch numbers, length, and diameter, calculated based on Horsfield ordering are slightly larger than Tawhai et al. [12] results. That is, the proposed model is not as symmetric as Tawhai’s model. In contrast to Strahler-based ratios, the Horsfield values of Rb, Rd, and Rl for the proposed model follow the cube-root relationship. The distribution of diameters of the dichotomous structure is related to the distribution of air flow in lung lobes based on a power law relationship [45]. The flow-dividing ratio (r ) plays a key role to compute the diameter proportional to the fluid flow. This parameter has been chosen randomly for each branch, but it has also been checked simultaneously whether the calculated diameter is consistent with the volume, which was supplied by that branch or not. Randomness of r makes the algorithm stochastic, that is, the generated trees are slightly different, but the geometrical and anatomical properties of all features of the trees are identical. In this study, we tried to calculate diameters and lengths independently through various generations. The diameters have been calculated according to the diameter-flow relationship mentioned before, and the lengths have been calculated based on the distance from the boundary. Kitaoka et al. [11] calculated lengths and diameters as interdependent parameters. They assigned a value to the length of each daughter branch that was three times of its diameter. Tawhai et al. [12] calculated the length of the branches of the whole structure based on the center of mass of the seed points which had been distributed randomly in host volume; then the diameters were randomly assigned to branches in each Horsfield order with arbitrarily chosen factors of 0.1. Tawhai et al. [13] estimated the diameters based on the cube-root relationship between branching ratio (Rb) and mean diameter ratio (Rd) in Strahler ordering. In two previous studies [12,13], diameters of branches with the same order are identical, which is not consistent with anatomic concepts. Whereas in the proposed method, the estimation of branch length is obtained according to the distance between each branch and its closest obstacle, including lobe boundary or other branches, in the path of the growing of the branch (Eq. 5). The advantage of this approach is that it can simulate the growth of the bronchial tree due to the enlargement of lung parenchyma during different ages. Although the length-to-diameter parameter is larger than published measurements [4,11,13], it can be fixed by constraining Γ max in different generation numbers based on anatomic data that are not currently available. Geometry analysis of the whole structure shows the branching pattern in the proposed method is slightly more asymmetric than in the previous presented model [6,12]. To generate a more symmetric branching pattern, the range of the variation of the flow-dividing ratio should be decreased and gets closer to the upper limit. Unfortunately, this would cause the branching

computer methods and programs in biomedicine 132 (2016) 1–10

process to be terminated faster, and the resulting model cannot supply the lung volume appropriately. It also reduces the branching angle. Variants of the value of the diameter exponent (n) according to the flow regime would be proposed as the other solution. However, this imposes the cost of a much more complex computation to the algorithm.

6.

Conclusion

We have shown a 3D model of the human conducting airway specific to an individual by the proposed algorithm that incorporates branching characteristics and space division. The structure has been divided into two segments, and the parametric OL-system was used to model each segment. Totally, 11 rewriting rules were provided to construct the bronchial tree. Each alphabetic character in rewriting rules played the role of a branch with its dimensional properties. The algorithm yields a spatially distributed model that is generally consistent with morphometric characteristics obtained by airway cast studies. The branching angle and the diameter are calculated by the equation relating diameter and air flow rate, which was presented by Kitaoka et al. [50]. In this study, the structure is generated separately into the lobes with accurate geometry, taken from the CT data. Furthermore, the lengths of branches are estimated based on its distance from the closest obstacle along the growth path (either lobe wall or other branches). This approach can be useful for simulating airways in lungs with structural abnormalities or simulating dynamical growth of the bronchial tree. Integrating a proper respiratory airway model and terminal branches will lead to a model which can be used to investigate the distribution of airflow in human respiratory system to design efficient drug delivery systems via inhalation.

Appendix. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.cmpb.2016.04.021.

REFERENCES

[1] M.G. Levitzky, Pulmonary Physiology, seventh ed., McGrawHill Medical, New York, 2007. [2] J.F. Murray, The Normal Lung: The Basis for Diagnosis and Treatment of Pulmonary Disease, second ed., Saunders, Philadelphia, 1986. [3] R.B. Schlesinger, L.A. McFadden, Comparative morphometry of the upper bronchial tree in six mammalian species, Anat. Rec. 199 (1981) 99–108. [4] E.R. Weibel, Morphometry of the Human Lung, Academic Press, New York, 1963. [5] K. Horsfield, G. Cumming, Morphology of the bronchial tree in man, J. Appl. Physiol. 24 (1968) 373–383. [6] K. Horsfield, G. Dart, D.E. Olson, G.F. Filley, G. Cumming, Models of the human bronchial tree, J. Appl. Physiol. 31 (1971) 207–217.

9

[7] K. Horsfield, A. Thurlbeck, Computer simulation of the geometry of the human bronchial tree, Bull. Math. Biol. 46 (1984) 389–398. [8] T.B. Martonen, Y. Yang, M. Dolovich, Definition of airway composition within gamma camera images, J. Thorac. Imaging 9 (1994) 188–197. [9] T.R. Nelson, D.K. Manchester, Modeling of lung morphogenesis using fractal geometries, IEEE Trans. Med. Imaging 7 (1988) 321–327. [10] R.F. Phalen, H.C. Yeh, G.M. Schum, O.G. Raabe, Application of an idealized model to morphometry of the mammalian tracheobronchial tree, Anat. Rec. 190 (1978) 167–176. [11] H. Kitaoka, R. Takaki, B. Suki, A three-dimensional model of the human airway tree, J. Appl. Physiol. 87 (1999) 2207– 2217. [12] M.H. Tawhai, A.J. Pullan, P.J. Hunter, Generation of an anatomically based three-dimensional model of the conducting airways, Ann. Biomed. Eng. 28 (2000) 793–802. [13] M.H. Tawhai, P. Hunter, J. Tschirren, J. Reinhardt, G. McLennan, E.A. Hoffman, CT-based geometry analysis and finite element models of the human and ovine bronchial tree, J. Appl. Physiol. 97 (2004) 2310–2321. [14] W. Park, E.A. Hoffman, M. Sonka, Segmentation of intrathoracic airway trees: a fuzzy logic approach, IEEE Trans. Med. Imaging 17 (1998) 489–497. [15] T.-Y. Law, P. Heng, Automated extraction of bronchus from 3D CT images of lung based on genetic algorithm and 3D region growing, Proc. SPIE Int. Soc. Opt. Eng. 3979 (2000) 906– 916. [16] D. Aykac, E.A. Hoffman, G. McLennan, J.M. Reinhardt, Segmentation and analysis of the human airway tree from three-dimensional X-ray CT images, IEEE Trans. Med. Imaging 22 (2003) 940–950. [17] P. Jiantao, C. Fuhrman, W.F. Good, F.C. Sciurba, D. Gur, A differential geometric approach to automated segmentation of human airway tree, IEEE Trans. Med. Imaging 30 (2011) 266–278. [18] P. Lo, B. van Ginneken, J.M. Reinhardt, T. Yavarna, P.A. de Jong, B. Irving, et al., Extraction of airways from CT (EXACT’09), IEEE Trans. Med. Imaging 31 (2012) 2093–2107. [19] N. Nowak, P. Kakade, A. Annapragada, Computational fluid dynamics simulation of airflow and aerosol deposition in human lungs, Ann. Biomed. Eng. 31 (2003) 374–390. [20] B. Ma, K. Lutchen, CFD simulation of aerosol deposition in an anatomically based human large–medium airway model, Ann. Biomed. Eng. 37 (2009) 271–285. [21] Y. Yin, J. Choi, E.A. Hoffman, M.H. Tawhai, C.-L. Lin, A multiscale MDCT image-based breathing lung model with time-varying regional ventilation, J. Comput. Phys. 244 (2013) 168–192. [22] N. Tgavalekos, J.G. Venegas, B. Suki, K.R. Lutchen, Relation between structure, function, and imaging in a threedimensional model of the lung, Ann. Biomed. Eng. 31 (2003) 363–373. [23] T. Law, P. Heng, Automatic centerline extraction for 3d virtual bronchoscopy, Proc. MICCAI 1935 (2000) 786–793. [24] A.P. Kiraly, W.E. Higgins, G. McLennan, E.A. Hoffman, J.M. Reinhardt, Three-dimensional human airway segmentation methods for clinical virtual bronchoscopy, Acad. Radiol. 9 (2002) 1153–1168. [25] T. Schlathoelter, C. Lorenz, I.C. Carlsen, S. Renisch, T. Deschamps, Simultaneous segmentation and tree reconstruction of the airways for virtual bronchoscopy, Proc. SPIE Int. Soc. Opt. Eng. 4684 (2002) 103–113. [26] A.P. Kiraly, J.P. Helferty, E.A. Hoffman, G. McLennan, W.E. Higgins, Three-dimensional path planning for virtual bronchoscopy, IEEE Trans. Med. Imaging 23 (2004) 1365– 1379.

10

computer methods and programs in biomedicine 132 (2016) 1–10

[27] A. Lindenmayer, Mathematical models for cellular interactions in development I. Filaments with one-sided inputs, J. Theor. Biol. 18 (1968) 280–299. [28] A. Lindenmayer, Mathematical models for cellular interactions in development II. Simple and branching filaments with two-sided inputs, J. Theor. Biol. 18 (1968) 300– 315. [29] P. Prusinkiewicz, A. Lindenmayer, The Algorithmic Beauty of Plants, Springer-Verlag, New York, 1990. [30] S. Ruoxi, J. Jinyuan, M. Jaeger, Intelligent tree modeling based on L-system, IEEE Int. Conf. CAID & CD (2009) 1096– 1100. [31] Y. Livny, F. Yan, M. Olson, B. Chen, H. Zhang, J. El-Sana, Automatic reconstruction of tree skeletal structures from point clouds, ACM Trans. Graph. 29 (2010) 1–8. [32] P. Prusinkiewicz, J. Hanan, R. Meˇch, An L-system-based plant modeling language, Proc. AGTIVE 1779 (2000) 395– 410. [33] O. Petrenko, M. Sbert, O. Terraz, D. Ghazanfarpour, 3Gmap L-systems grammar application to the modeling of flowering plants, Intel. Comput. Graph. 441 (2013) 1– 21. [34] L. Xin, L. Xu, D. Li, D. Fu, The 3D reconstruction of greenhouse tomato plant based on real organ samples and parametric L-system, Proc. SPIE ICDIP 9159 (2014) 915904– 915905. [35] M. Zamir, Arterial branching within the confines of fractal l-system formalism, J. Gen. Physiol. 118 (2001) 267–276. [36] X. Liu, H. Liu, A. Hao, Q. Zhao, Simulation of blood vessels for surgery simulators, Int. Conf. MVHI (2010) 377–380. [37] A. Lindenmayer, Developmental systems without cellular interactions, their languages and grammars, J. Theor. Biol. 30 (1971) 455–484. [38] EXACT’09, Extraction of Airways from CT 2009. http://image .diku.dk/exact/, 2009.

[39] P. Nardeli, P. Cantillon-Murphy, Semi-automated airway segmentation for lung CT datasets, Proc. Int. Congr. Exhib. CARS (2013). [40] 3D Slicer, A multi-platform, free and open source software package for visualization and medical image computing. http://www.slicer.org/. [41] Y. Gao, A. Tannenbaum, R. Kikinis, Simultaneous multiobject segmentation using local robust statistics and contour interaction, Proc. MICCAI, Workshop on MCV 6533 (2011) 195–203. [42] C.D. Murray, The physiological principle of minimum work, Proc. Natl. Acad. Sci. U.S.A. 12 (1926) 207–214. [43] B.B. Mandelbrot (Ed.), The Fractal Geometry of Nature, W.H. Freeman, San Francisco, 1983. Updated and augmented. [44] H.B.M. Uylings, Optimization of diameters and bifurcation angles in lung and vascular tree structures, Bull. Math. Biol. 39 (1977) 509–520. [45] H. Kitaoka, B. Suki, Branching design of the bronchial tree based on a diameter-flow relationship, J. Appl. Physiol. 82 (1997) 968–976. [46] O.G. Raabe, H.-C. Yeh, G.M. Schum, R.F. Phalen, Tracheobronchial Geometry: Human, Dog, Rat, Hamster, 1976. [47] B. Sapoval, M. Filoche, E.R. Weibel, Smaller is better—but not too small: a physical scale for the design of the mammalian pulmonary acinus, Proc. Natl. Acad. Sci. U.S.A. 99 (2002) 10411–10416. [48] K. Horsfield, G. Cumming, Angles of branching and diameters of branches in the human bronchial tree, Bull. Math. Biophys. 29 (1967) 245–259. [49] B. Haefeli-Bleuer, E.R. Weibel, Morphometry of the human pulmonary acinus, Anat. Rec. 220 (1988) 401–414. [50] K. Horsfield, F.G. Relea, G. Gumming, Diameter, length and branching ratios in the bronchial tree, Respir. Physiol. 26 (1976) 351–356.