TY - JOUR
AU - Moustakas, Konstantinos
AB - 1 Introduction In the past years, a multitude of studies paves the way for the generation of patient-specific computational models of lung structure and function. Early studies focused on airway morphometry generating the first human bronchial trees models [1]. These studies employed casts to decipher the relationship between bronchi lengths, branching angles and airway diameters [2]. On this basis, researchers built and validated a simulation model of airway morphogenesis from generation 1 to generation 23 [3, 4]. Deterministic parameterized bronchial tree generation algorithms used as single input the location of the first one or two generations and the lung volume, extracted directly from computed tomography, thus constituting the core of patient-specific modelling [5–7]. Personalized boundary conditions based on diagnostic imaging were combined with generative approaches and lumped models of resistive trees [8, 9] constituting the state of the art in pulmonary system modelling. Later studies incorporated patient-specific boundary conditions into computational fluid dynamics to examine flow regimes, wall stresses and aerosol deposition. In the same direction, modelling the airflow in cases of constrictive conditions, such as asthma and chronic obstructive pulmonary disease (COPD) became feasible with the aforementioned approaches. Wall constriction and remodelling combined with patient-specific boundary conditions allowed the quantification of breathing conditions for asthmatic patients. Motivated by these advancements we introduce an end-to-end modelling approach that produces Anatomically Valid Airway tree conformations(AVATREE). Such conformations are adapted to personalized geometry and boundary conditions derived from diagnostic imaging and well-established airway extraction methods. Specifically, this study aims to provide an open-source simulation framework to (i) exploit imaging data so as to provide patient-specific representations (ii) perform structural analysis (iii) extend the segmented airway tree to predict the airway branching across the whole lung volume (iv) visualize probabilistic confidence maps of generation data (v) simulate bronchoconstriction to (vi) access patient-specific airway functionality (vii) perform fluid dynamics simulation in patient-specific boundary conditions to access pulmonary function. 1.1 Background & related work While early studies focused mainly on quantitative modeling approaches to gain insight into the lung function without an explicit link to the lung’s structure, with the advancements in computing power and the current medical imaging capabilities, the interest in the simulation of lung function based on personalized geometric models that incorporate the essential structural features of the lungs, has significantly increased [10]. By now, many studies propose the development and adoption of mathematical and geometrical models to study the structure of the airways and pulmonary physiology. Some address the problem of airway tree segmentation from CT images, while others analyze the branching patterns and bifurcations through airway morphometry or mathematical modelling. In this section, we briefly present representative approaches and introduce definitions for the different computational steps required during airway and pulmonary structure modelling and simulation. Airway segmentation, bronchial morphometry and tree branching, mathematical models of bifurcating distributing systems are required to derive patient-specific structural and functional modelling approaches. Early studies on airway morphometry [11–13] used casts of human lungs to study branching patterns and the relation between airway lengths and diameters. The most commonly used conducting airway model has been Weibel’s symmetric model “A” [14]. The airway position has also been described by Horsfield order [1] and Strahler order [15]. Later on, with the advancement of medical imaging techniques, the extraction of airway structure and lung volume from imaging started to play an important role in the analysis of pulmonary diseases. A literature review on the analysis of lung CTs, including segmentation of the various pulmonary structures, can be found in [16], while a comparative study of automated and semi-automated segmentation methods of the airway tree from CT images was presented in [17]. Overall, segmentation approaches can be classified into methods based on morphology [18], morphological aggregation [19], voxel classification [20], adaptive region growing with constraints [21–28], tube similarity [29, 30] and gradient vector flow [31]. Several implementations of the aforementioned approaches are available in the literature. The tube segmentation framework [30] utilizing gradient vector flow [31] and the FAST heterogeneous medical image computing and visualization framework [32] utilizing the seeded region growing approach. AVATREE employs airway segmentation as a first step to obtain the personalized structure in the first generations, while the more advanced generations are simulated based on a tree extension algorithm. Furthermore, mathematical models of the airway structure were formulated to derive branching and structural rules. Deterministic mathematical models of bifurcating distributing systems were examined [3] setting the basis for modelling bronchial tree branching as a function of available lung space [4]. Deriving airway diameter as a relation of branching features facilitates full determination of the geometry given skeletal representations. Several studies mention scaling properties [33–36] for the airway diameters so that the average diameter D of a given airway at generation G is the product of the diameter of the trachea D0. Furthermore, Kamiya et al. [37] validated the relation between airway diameter and branching angles and, Kitaoka et al. [2] proposed a branching model allowing the prediction of the relationship between branching angle and flow rate and between airway length and diameter. Experimental studies verified the validity of the aforementioned methods. For the surface reconstruction of airway surface Tawhai et al. [5, 6, 38] employed fitting cubic Hermite surfaces as described in [39]. The study of Hegedus et al. [40] generated surface models of idealized bifurcation through mathematical modelling rigorously extending the previous definitions [41]. The aforementioned studies are relevant to our approach. To avoid the definition of special rules in the reconstruction of the surface of bifurcations we define the same boundary conditions, as the Poisson-reconstructed surface of a sampled point cloud. Towards patient-specific structural and functional modelling, Tawhai et al. and Lin et al. [6, 42] studied the imposition of patient-specific boundary conditions to generate 1-dimensional and three-dimensional computational models taking into consideration the effects of turbulence. Towards the same direction, a review article [10] provides insight into multiscale finite element models of lung structure and function aiming towards a computational framework for bridging the spatial scales from molecular to the whole organ. Bordas et al. [7] developed an image analysis and modelling pipeline applied to healthy and asthmatic patient scans to produce complete personalized airway models to the acinar level incorporating CT acquisition, lung and lobar segmentation, airway segmentation and centerline extraction, algorithmic generation of distal airways and zero-dimensional models. Their implementation and results were included into Chaste framework [43], an open-source framework to facilitate computational modelling in heart, lung and soft tissue simulations. Towards the same direction, Montesantos et al. [44] presented a detailed algorithm for the generation of an individualized 3D deterministic model of the conducting part of the human tracheobronchial tree. With respect to the aforementioned studies, our work focuses on generating surface meshes of extended patient-based bronchial trees, suitable for computational fluid dynamics (CFD) simulations, along with a toolbox to simulate constriction of the airways. Several authors employed CFD to investigate flow regimes in the human lung. In our previous work [45, 46] we performed narrowing deformations in CT extracted lung geometries to simulate constrictive conditions. Other studies in the same category include simulations for CT-based patient specific geometries [47–52], particle deposition [53–56], constrictive pulmonary diseases [45, 46, 57–59], micro-airway flow regimes [57], turbulence modelling [60], four-dimensional (space and time) dynamic simulations [61], ventilation heterogeneity [57], airflow in the acinar region [62]. Validation studies conducted by Montesantos et al. [63] include morphometric studies on healthy and asthmatic patients providing among others, measurements of branching angles, length and diameter of airways as a function of generation. Such measurements are employed by our study for macroscopic validation of the generated trees. 1.2 Motivation and contributions The objective in this field of research is to enable the prediction of gas flow [51, 55], gas mixing [64], heat transfer [65], particle deposition [46, 54, 66, 67], and ventilation distribution [68] in the pulmonary system. Lung ventilation patterns prediction [69, 70] can provide grounds for performance and fatigue estimation in high-frequency ventilation cases [71], disease severity quantification, such as in asthma and COPD, and give insight into drug delivery or even in transfer of harmful particulates. Motivated by the recent advances in this field and building upon previous work [46], we developed an end-to-end approach facilitating pulmonary structural modelling that is based on the definition of the personalised boundary conditions required for fluid dynamics simulations. Specifically, in this work we present an open-source simulation framework that utilizes imaging data to provide patient-based representations of the structural models of the bronchial tree, build and extend 1-dimensional graph representations of the bronchial tree, generate 3D surface models of extended bronchial tree models appropriate for CFD simulations generate probabilistic visualization of airway generations projected on the personalized CT imaging data, 2 perform validation studies and provide comparison with relevant state-of-the-art approaches provide an open-source toolbox in C++ and a graphical user interface integrating modelling functionalities. The rest of the paper is organized as follows. Section 2 analyzes the individual components of AVATREE, Section 3 commends on the results of our approach while Section 4 concludes this paper. 1.1 Background & related work While early studies focused mainly on quantitative modeling approaches to gain insight into the lung function without an explicit link to the lung’s structure, with the advancements in computing power and the current medical imaging capabilities, the interest in the simulation of lung function based on personalized geometric models that incorporate the essential structural features of the lungs, has significantly increased [10]. By now, many studies propose the development and adoption of mathematical and geometrical models to study the structure of the airways and pulmonary physiology. Some address the problem of airway tree segmentation from CT images, while others analyze the branching patterns and bifurcations through airway morphometry or mathematical modelling. In this section, we briefly present representative approaches and introduce definitions for the different computational steps required during airway and pulmonary structure modelling and simulation. Airway segmentation, bronchial morphometry and tree branching, mathematical models of bifurcating distributing systems are required to derive patient-specific structural and functional modelling approaches. Early studies on airway morphometry [11–13] used casts of human lungs to study branching patterns and the relation between airway lengths and diameters. The most commonly used conducting airway model has been Weibel’s symmetric model “A” [14]. The airway position has also been described by Horsfield order [1] and Strahler order [15]. Later on, with the advancement of medical imaging techniques, the extraction of airway structure and lung volume from imaging started to play an important role in the analysis of pulmonary diseases. A literature review on the analysis of lung CTs, including segmentation of the various pulmonary structures, can be found in [16], while a comparative study of automated and semi-automated segmentation methods of the airway tree from CT images was presented in [17]. Overall, segmentation approaches can be classified into methods based on morphology [18], morphological aggregation [19], voxel classification [20], adaptive region growing with constraints [21–28], tube similarity [29, 30] and gradient vector flow [31]. Several implementations of the aforementioned approaches are available in the literature. The tube segmentation framework [30] utilizing gradient vector flow [31] and the FAST heterogeneous medical image computing and visualization framework [32] utilizing the seeded region growing approach. AVATREE employs airway segmentation as a first step to obtain the personalized structure in the first generations, while the more advanced generations are simulated based on a tree extension algorithm. Furthermore, mathematical models of the airway structure were formulated to derive branching and structural rules. Deterministic mathematical models of bifurcating distributing systems were examined [3] setting the basis for modelling bronchial tree branching as a function of available lung space [4]. Deriving airway diameter as a relation of branching features facilitates full determination of the geometry given skeletal representations. Several studies mention scaling properties [33–36] for the airway diameters so that the average diameter D of a given airway at generation G is the product of the diameter of the trachea D0. Furthermore, Kamiya et al. [37] validated the relation between airway diameter and branching angles and, Kitaoka et al. [2] proposed a branching model allowing the prediction of the relationship between branching angle and flow rate and between airway length and diameter. Experimental studies verified the validity of the aforementioned methods. For the surface reconstruction of airway surface Tawhai et al. [5, 6, 38] employed fitting cubic Hermite surfaces as described in [39]. The study of Hegedus et al. [40] generated surface models of idealized bifurcation through mathematical modelling rigorously extending the previous definitions [41]. The aforementioned studies are relevant to our approach. To avoid the definition of special rules in the reconstruction of the surface of bifurcations we define the same boundary conditions, as the Poisson-reconstructed surface of a sampled point cloud. Towards patient-specific structural and functional modelling, Tawhai et al. and Lin et al. [6, 42] studied the imposition of patient-specific boundary conditions to generate 1-dimensional and three-dimensional computational models taking into consideration the effects of turbulence. Towards the same direction, a review article [10] provides insight into multiscale finite element models of lung structure and function aiming towards a computational framework for bridging the spatial scales from molecular to the whole organ. Bordas et al. [7] developed an image analysis and modelling pipeline applied to healthy and asthmatic patient scans to produce complete personalized airway models to the acinar level incorporating CT acquisition, lung and lobar segmentation, airway segmentation and centerline extraction, algorithmic generation of distal airways and zero-dimensional models. Their implementation and results were included into Chaste framework [43], an open-source framework to facilitate computational modelling in heart, lung and soft tissue simulations. Towards the same direction, Montesantos et al. [44] presented a detailed algorithm for the generation of an individualized 3D deterministic model of the conducting part of the human tracheobronchial tree. With respect to the aforementioned studies, our work focuses on generating surface meshes of extended patient-based bronchial trees, suitable for computational fluid dynamics (CFD) simulations, along with a toolbox to simulate constriction of the airways. Several authors employed CFD to investigate flow regimes in the human lung. In our previous work [45, 46] we performed narrowing deformations in CT extracted lung geometries to simulate constrictive conditions. Other studies in the same category include simulations for CT-based patient specific geometries [47–52], particle deposition [53–56], constrictive pulmonary diseases [45, 46, 57–59], micro-airway flow regimes [57], turbulence modelling [60], four-dimensional (space and time) dynamic simulations [61], ventilation heterogeneity [57], airflow in the acinar region [62]. Validation studies conducted by Montesantos et al. [63] include morphometric studies on healthy and asthmatic patients providing among others, measurements of branching angles, length and diameter of airways as a function of generation. Such measurements are employed by our study for macroscopic validation of the generated trees. 1.2 Motivation and contributions The objective in this field of research is to enable the prediction of gas flow [51, 55], gas mixing [64], heat transfer [65], particle deposition [46, 54, 66, 67], and ventilation distribution [68] in the pulmonary system. Lung ventilation patterns prediction [69, 70] can provide grounds for performance and fatigue estimation in high-frequency ventilation cases [71], disease severity quantification, such as in asthma and COPD, and give insight into drug delivery or even in transfer of harmful particulates. Motivated by the recent advances in this field and building upon previous work [46], we developed an end-to-end approach facilitating pulmonary structural modelling that is based on the definition of the personalised boundary conditions required for fluid dynamics simulations. Specifically, in this work we present an open-source simulation framework that utilizes imaging data to provide patient-based representations of the structural models of the bronchial tree, build and extend 1-dimensional graph representations of the bronchial tree, generate 3D surface models of extended bronchial tree models appropriate for CFD simulations generate probabilistic visualization of airway generations projected on the personalized CT imaging data, 2 perform validation studies and provide comparison with relevant state-of-the-art approaches provide an open-source toolbox in C++ and a graphical user interface integrating modelling functionalities. The rest of the paper is organized as follows. Section 2 analyzes the individual components of AVATREE, Section 3 commends on the results of our approach while Section 4 concludes this paper. 2 Materials and methods The processing pipeline uses as input CT images and is presented in Fig 1. Airway segmentation is applied to the imaging data to extract a 3D surface mesh and a 1-dimensional representation of the airways. We employ the extended 1D graph to derive visualization of probabilistic airway generation labels in the space of the subjects’ anatomy as defined by the CT images and to generate a 3D surface defining personalized boundary conditions, that can be employed as input for computational fluid dynamics simulations. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 1. Processing pipeline of AVATREE. https://doi.org/10.1371/journal.pone.0230259.g001 2.1 Segmentation and airways centerline extraction The input of the presented approach is unlabeled CT scans required to extract bronchial tree and airways structural features. For the definition of the lung volume, CT-based lung segmentation and annotation is required. For lung segmentation we employ the FAST heterogeneous medical image computing and visualization framework [32] is employed. The result of lung segmentation process is a binary mask visualized in Fig 1. As a next step, we perform further processing of the segmentation result to distinguish left and right lungs. The process is described below: A second region growing takes place starting from a single random point inside any of the segmented region only if all its immediate neighbours bare the same label. To advance the region growing front, all points neighbouring a candidate voxel must not include background voxel. This region is given a new label. Steps 1 and 2 are repeated for the other lung volume. The result is an image with three labels(background and two lung volumes). To distinguish left or right we employ the directed graph extracted from the main airways and follow the generic rule according to which the topological distance the topological distance between the bifurcations of the first and the second generation is longer in the left lung. The next step involves the segmentation of the first generations of the airways that are identifiable in the patient’s CT image, but any available airway tree segmentation method can be also applied. For this purpose we investigated two algorithms. The first algorithm is the gradient vector flow [29, 31] which achieved high accuracy with low false-positive rate (only 1.44%) in a comparative study [17] in the context of the EXACT09 airway segmentation challenge. The second is a standard and stable approach based on seeded region growing [72]. The former is included in the tube segmentation framework [30] and the latter in FAST heterogeneous medical image computing and visualization framework [32]. Let’s denote with I(x), , the gray level 3D medical image, where x = (x, y, z), x ∈ Ω is a voxel in the spatial domain of the volumetric imaging data. The output of the segmentation algorithm for the airways is a binary image SA of equal size with I. Likewise, the output of the segmentation algorithm for the lung volumes is a binary image SL of equal size with I. The result is presented in Figs 1 and 2, and utilized to generate prediction of full bronchial tree structures based on personalized lung volumes. To derive the centerline from SA a multitude of methods is provided in the literature including skeletonization or thinning. Fig 2 presents the up-to-four generations centerline of the airways. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 2. Extraction of airway surface and centerline. https://doi.org/10.1371/journal.pone.0230259.g002 This 1D representation of the bronchial tree is modelled by an undirected graph where is the set of vertices and is the set of edges. Each vertex, indexed by i, can be represented as a point vi = (xi, yi, zi). We denote the function N(vi) yielding the set of vertices indices neighbouring vertex i. The undirected graph is extracted by FAST framework [32] and converted into a directed graph with the following process. Initially, the graph starting point is defined as the one closest to the air inlet, i.e. the oral cavity or the trachea. Given index y the starting point for we generate the directed tree . We define as distal point the vertex of the graph with no children and distal branch the edge containing a distal point. Algorithm 1 Detection of inlet in undirected graph Input: Graph Output: Index of graph inlet y 1 procedure Derivation of graph inlet 2  Initialize set 3  for each vertex vi do 4   if |N(vi)| > 2 then 5    for each n ∈ N(vi) do 6     Initialize empty set 7     while N(vn)<3 do 8      for each m ∈ N(vn) do 9       if then 10      11   12   2.2 Generation of extended bronchial tree Since higher generations cannot be identified from the personal imaging data, we extend the bronchial tree based on population-wise empirical observations. Initially the directed graph generated by the procedure explained in subsection 2.1 is pruned. Specifically, the extracted 1-dimensional representation is processed to include all the bifurcations located at the end of a given generation so as to facilitate the volume filling algorithm. Fig 1 shows the result of pruning where all generations after the nth have been pruned. The corrected tree is subsequently used for the bronchial tree extension. The generation process utilizes the bronchial tree extension algorithm initially proposed by Tawhai et al. [4] and later improved by Bordas et al. [7] while introducing a few safeguards to allow maximal space utilization. The bronchial tree extension algorithm can be described by the following steps. For each lung subvolume and : Generate a point cloud sampling the subvolume with a uniform random process. Fig 3 depicts the uniform sampling of each lung subvolume with a total number of n = 30000 points [4, 6]. Assign a seed point to the closest distal branch as presented in Fig 3. Calculate the center of mass of the sampled points as presented in Fig 3. (1) Employ principal component analysis (PCA) on the set of sampled points to define the splitting plane. The motivation for employing PCA is to address a space utilization aspect. The direction of the eigenvector with the greatest norm indicates the dimension of the data with the greatest variance denoting the direction where more space is available for the branches to grow. Picking a plane so that the resulting bounding box demonstrates the lowest possible variation, inhibits the appearance of very long branches. Given data points D = [p1 p1 p1⋯pn], A = DDT is the auto-correlation matrix. Direct singular value decomposition yields A = UΣUT where U = [u1 u2 u3]. Then the largest eigenvector is defined as um = max1≤i≤3 ui. Given the vector d expressing the direction of the distal airway, the splitting plane is described by center of mass c and vector d × [d × um]. The selected plane maximizes the available space for each new subdivision. Fig 4 presents a splitting plane splitting the set of points into two subdivisions. Calculate the centroid of each new subdivision. For each centroid define line segment starting from seed point extending 40% of the distance towards centroid of the subdivision. If a newly created branch is smaller that 2mm, it is considered as terminal. The process is repeated until no seed points remain. Any branch found outside the lung volume is removed along with children branches. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 3. Definition of center of mass (orange dot) and distal branches (yellow dots). https://doi.org/10.1371/journal.pone.0230259.g003 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 4. Definition of splitting plane for bifurcating distributive structures. https://doi.org/10.1371/journal.pone.0230259.g004 It is important to denote that the presented pipeline enables the generation of a tunable user-defined number of generations. If n is the number of desired generations, we set stopping criteria, in the extension of the bronchial tree until 2(n+1) bifurcations have been reached. The resulting 1D representation (Fig 5) predicts the location of the bifurcating distributive [3] structure given the patient-specific available space. The outcome of the volume filling algorithm will be used later to create maps that express the probability of a voxel to belong to a certain generation. This information when projected on CT slices can be a very informative and powerful clinical decision support tool. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 5. Extended bronchial tree (a) 9 generations, (b) 12 generations. https://doi.org/10.1371/journal.pone.0230259.g005 2.3 Spatial probability maps of branching properties The location of each new generation branch is calculated as explained before and provides a random sample out of all possible bronchial tree conformations. In this step of the proposed framework we produce probabilistic maps for each generation branch that provide estimates of the spatial probability to encounter a certain generation at some point of the imaging data. Such a probabilistic model allows to optimize clinical decision making by accounting for the branches’ distributional uncertainty. Let’s denote with Wg(x), , the probability map for generation g where x = (x, y, z), x ∈ Ω is a voxel in the spatial domain of the volumetric imaging data. Then (2) where d is the distance of voxel x to the closest edge of labeled with generation g. Parameter σ is set experimentally to σ = 1. The extracted spatial map is overlaid on the CT scans, as shown in Fig 6, providing insightful visualization of spatial likelihood for each branching generation. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 6. Visualization of fourth generation probability map. First column corresponds to axial view. Second column corresponds to coronal view. Third column corresponds to sagital view. The first row depicts raw imaging data, the second row presents the probabilistic maps for second generation, the third row presents the probabilistic maps for second generation and the fourth row presents the probabilistic maps for second generation. https://doi.org/10.1371/journal.pone.0230259.g006 2.4 Surface generation on predicted 1-dimensional representation The extracted 3D geometries are required to conduct studies on computational fluid dynamics, particle transfer and deposition, ventilation, stress analysis and deformation simulations. Marching cubes algorithm [73] is a very well established method implemented in FAST [32] allowing the generation of 3D geometric models from airway segmentation label maps. The constriction simulation method aims to generate 3D tubular surface structures with smaller diameters. To this end, Laplacian surface contraction offers a solution that deforms the geometry pushing the vertices towards the direction of the inward normals. The extension of the extracted centerline generates a predictive representation of the bronchial tree given the available space. However, for the outcome of space-filling algorithms to be useful in fluid dynamics simulation, particle deposition simulation or stress finite element analysis based studies, the boundary conditions in the form of triangular 3D meshes need to be defined. Initially, as a simplified approach, to define the diameter of each generation we can employ the power law consistent with Murray’s law of symmetric branching [33, 34]. (3) where d0 denotes the branch diameter of the trachea and dz the branch diameter for generation z. Furthermore, if we take into account that each branch demonstrates different branching angle and diameter properties, the relation between airway diameter (d) and branching angles (θ) is based on the following rules validated by Kamiya et al. [37] and Kitaoka et al. [2]: (4) (5) where the index 0 stands for the parent branch, and the indices 1 and 2 for the two children branches, respectively. To reconstruct the lung surface we employ a point cloud sampling approach as input for Poisson surface reconstruction. The outcome is a smooth surface with smooth transitions instead of abrupt transitions in the intersection with the original tubular meshes. The tubular-shaped point cloud is sampled using a uniform random distribution. A clean-up step, visualized in Fig 7, ensures that no point can be found in distance less than the prescribed diameter of every available branch. The resulting point cloud is used to compute normals. A bilateral normal smoothing [74] function prepares the point cloud for Poisson surface reconstruction [75]. smoothing the point normals. This step facilitates the surface reconstruction in bifurcations and transitional parts. Furthermore, since the directed graph is extracted where each point on the centerline corresponds to a point on the lung surface it is possible to further deform the surface with a custom function or pattern. The generated surface for seven and ten generations is presented in Fig 8. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 7. Outcome of point cloud generation process. Extra refinements remove the inner points facilitating normal estimation for Poisson surface reconstruction. https://doi.org/10.1371/journal.pone.0230259.g007 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 8. Generation of surface from the extended bronchial tree centerline. https://doi.org/10.1371/journal.pone.0230259.g008 2.5 Simulation of constrictive pulmonary diseases’ effect on airway tree This section aims to provide the methodology for simulation of broncho-constriction allowing to subsequently estimate the dynamic behaviour of the lung airways in the case of an exacerbation event. A bronchial tree 3D geometry is the input for this process yielding as output contracted airways. The proposed geometry contraction procedure is presented by Nousias et al. [45] and Lalas et al. [46] and is an extension of the work of Au et al. [76] facilitating a Laplacian smoothing process that shifts vertices along the estimated curvature normal direction. The airway geometric model consists of connected triangles forming the boundary conditions. Each triangular mesh can be described as where is the set of vertices, is the set of edges and is the set of faces constituting the 3D surface. Each vertex i can be represented as a point vi = (xi, yi, zi), ∀i = 1, 2, ⋯, N. For each face fi, ∀i = 1, 2, ⋯, l we denote the centroid (6) The outward unit normal to the face fi (located at the centroid mi) is calculated as : (7) where are the vertices corresponding to face fi. Given the curvature flow Laplacian operator, the product δ = LV approximates the inward curvature flow normals [77]. The motivation for employing the curvature flow Laplacian operator [78] on the mesh is that its output is not affected by mesh density. Specifically, (8) where Ai is the one-ring area, κi is the local curvature and ni is the inward curvature flow normal of the ith vertex. The positions of the vertices satisfying LV = 0 result in a zero volume string-like mesh and can be used to simulate mesh contraction. However, since such an optimization problem has more than one solutions, further constraints are required [76]. Introducing weighting matrices can smoothly move vertex positions towards the direction of the inward unit normal by iteratively solving the following minimization problem (9) where corresponds to the vertex positions before the contraction at each iteration. The weighting matrices WH and WL regulate the mesh contraction and mesh attraction, respectively. Initially, we set them to and WH = I, where is the identity matrix, k a double constant experimentally tuned to 10−3 and A the average face area of the model. Eq (9) can be expressed as (10) The analytical solution can be formulated as (11) where matrices A and b are given by (12) After each iteration t we update the contraction and inflation weights to be used in iteration t + 1 so that and , where and are the original and the current one-ring area respectively. The Laplacian matrix for iteration t + 1, Lt+1 is also recomputed. On these grounds, to simulate broncho-constriction we require to reduce the airway diameter to a predefined level of narrowing is reached. This level is defined by certain termination criteria [45, 46]. Thus, a metric is required that measures the diameter of the bronchi under process. To estimate the degree of contraction of the airway’s geometry after each iteration, we employ a shape diameter function (SDF) based scheme [79] implemented in [80] that evaluates the local volume based on the estimated local diameter assigned to each face of the mesh, also known as raw SDF values. Measuring the volume before and after the Laplacian contraction iteration can set the termination criteria. Fig 9 presents a simulation of constrictive pulmonary conditions. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 9. Simulation of constrictive pulmonary conditions. https://doi.org/10.1371/journal.pone.0230259.g009 2.1 Segmentation and airways centerline extraction The input of the presented approach is unlabeled CT scans required to extract bronchial tree and airways structural features. For the definition of the lung volume, CT-based lung segmentation and annotation is required. For lung segmentation we employ the FAST heterogeneous medical image computing and visualization framework [32] is employed. The result of lung segmentation process is a binary mask visualized in Fig 1. As a next step, we perform further processing of the segmentation result to distinguish left and right lungs. The process is described below: A second region growing takes place starting from a single random point inside any of the segmented region only if all its immediate neighbours bare the same label. To advance the region growing front, all points neighbouring a candidate voxel must not include background voxel. This region is given a new label. Steps 1 and 2 are repeated for the other lung volume. The result is an image with three labels(background and two lung volumes). To distinguish left or right we employ the directed graph extracted from the main airways and follow the generic rule according to which the topological distance the topological distance between the bifurcations of the first and the second generation is longer in the left lung. The next step involves the segmentation of the first generations of the airways that are identifiable in the patient’s CT image, but any available airway tree segmentation method can be also applied. For this purpose we investigated two algorithms. The first algorithm is the gradient vector flow [29, 31] which achieved high accuracy with low false-positive rate (only 1.44%) in a comparative study [17] in the context of the EXACT09 airway segmentation challenge. The second is a standard and stable approach based on seeded region growing [72]. The former is included in the tube segmentation framework [30] and the latter in FAST heterogeneous medical image computing and visualization framework [32]. Let’s denote with I(x), , the gray level 3D medical image, where x = (x, y, z), x ∈ Ω is a voxel in the spatial domain of the volumetric imaging data. The output of the segmentation algorithm for the airways is a binary image SA of equal size with I. Likewise, the output of the segmentation algorithm for the lung volumes is a binary image SL of equal size with I. The result is presented in Figs 1 and 2, and utilized to generate prediction of full bronchial tree structures based on personalized lung volumes. To derive the centerline from SA a multitude of methods is provided in the literature including skeletonization or thinning. Fig 2 presents the up-to-four generations centerline of the airways. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 2. Extraction of airway surface and centerline. https://doi.org/10.1371/journal.pone.0230259.g002 This 1D representation of the bronchial tree is modelled by an undirected graph where is the set of vertices and is the set of edges. Each vertex, indexed by i, can be represented as a point vi = (xi, yi, zi). We denote the function N(vi) yielding the set of vertices indices neighbouring vertex i. The undirected graph is extracted by FAST framework [32] and converted into a directed graph with the following process. Initially, the graph starting point is defined as the one closest to the air inlet, i.e. the oral cavity or the trachea. Given index y the starting point for we generate the directed tree . We define as distal point the vertex of the graph with no children and distal branch the edge containing a distal point. Algorithm 1 Detection of inlet in undirected graph Input: Graph Output: Index of graph inlet y 1 procedure Derivation of graph inlet 2  Initialize set 3  for each vertex vi do 4   if |N(vi)| > 2 then 5    for each n ∈ N(vi) do 6     Initialize empty set 7     while N(vn)<3 do 8      for each m ∈ N(vn) do 9       if then 10      11   12   2.2 Generation of extended bronchial tree Since higher generations cannot be identified from the personal imaging data, we extend the bronchial tree based on population-wise empirical observations. Initially the directed graph generated by the procedure explained in subsection 2.1 is pruned. Specifically, the extracted 1-dimensional representation is processed to include all the bifurcations located at the end of a given generation so as to facilitate the volume filling algorithm. Fig 1 shows the result of pruning where all generations after the nth have been pruned. The corrected tree is subsequently used for the bronchial tree extension. The generation process utilizes the bronchial tree extension algorithm initially proposed by Tawhai et al. [4] and later improved by Bordas et al. [7] while introducing a few safeguards to allow maximal space utilization. The bronchial tree extension algorithm can be described by the following steps. For each lung subvolume and : Generate a point cloud sampling the subvolume with a uniform random process. Fig 3 depicts the uniform sampling of each lung subvolume with a total number of n = 30000 points [4, 6]. Assign a seed point to the closest distal branch as presented in Fig 3. Calculate the center of mass of the sampled points as presented in Fig 3. (1) Employ principal component analysis (PCA) on the set of sampled points to define the splitting plane. The motivation for employing PCA is to address a space utilization aspect. The direction of the eigenvector with the greatest norm indicates the dimension of the data with the greatest variance denoting the direction where more space is available for the branches to grow. Picking a plane so that the resulting bounding box demonstrates the lowest possible variation, inhibits the appearance of very long branches. Given data points D = [p1 p1 p1⋯pn], A = DDT is the auto-correlation matrix. Direct singular value decomposition yields A = UΣUT where U = [u1 u2 u3]. Then the largest eigenvector is defined as um = max1≤i≤3 ui. Given the vector d expressing the direction of the distal airway, the splitting plane is described by center of mass c and vector d × [d × um]. The selected plane maximizes the available space for each new subdivision. Fig 4 presents a splitting plane splitting the set of points into two subdivisions. Calculate the centroid of each new subdivision. For each centroid define line segment starting from seed point extending 40% of the distance towards centroid of the subdivision. If a newly created branch is smaller that 2mm, it is considered as terminal. The process is repeated until no seed points remain. Any branch found outside the lung volume is removed along with children branches. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 3. Definition of center of mass (orange dot) and distal branches (yellow dots). https://doi.org/10.1371/journal.pone.0230259.g003 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 4. Definition of splitting plane for bifurcating distributive structures. https://doi.org/10.1371/journal.pone.0230259.g004 It is important to denote that the presented pipeline enables the generation of a tunable user-defined number of generations. If n is the number of desired generations, we set stopping criteria, in the extension of the bronchial tree until 2(n+1) bifurcations have been reached. The resulting 1D representation (Fig 5) predicts the location of the bifurcating distributive [3] structure given the patient-specific available space. The outcome of the volume filling algorithm will be used later to create maps that express the probability of a voxel to belong to a certain generation. This information when projected on CT slices can be a very informative and powerful clinical decision support tool. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 5. Extended bronchial tree (a) 9 generations, (b) 12 generations. https://doi.org/10.1371/journal.pone.0230259.g005 2.3 Spatial probability maps of branching properties The location of each new generation branch is calculated as explained before and provides a random sample out of all possible bronchial tree conformations. In this step of the proposed framework we produce probabilistic maps for each generation branch that provide estimates of the spatial probability to encounter a certain generation at some point of the imaging data. Such a probabilistic model allows to optimize clinical decision making by accounting for the branches’ distributional uncertainty. Let’s denote with Wg(x), , the probability map for generation g where x = (x, y, z), x ∈ Ω is a voxel in the spatial domain of the volumetric imaging data. Then (2) where d is the distance of voxel x to the closest edge of labeled with generation g. Parameter σ is set experimentally to σ = 1. The extracted spatial map is overlaid on the CT scans, as shown in Fig 6, providing insightful visualization of spatial likelihood for each branching generation. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 6. Visualization of fourth generation probability map. First column corresponds to axial view. Second column corresponds to coronal view. Third column corresponds to sagital view. The first row depicts raw imaging data, the second row presents the probabilistic maps for second generation, the third row presents the probabilistic maps for second generation and the fourth row presents the probabilistic maps for second generation. https://doi.org/10.1371/journal.pone.0230259.g006 2.4 Surface generation on predicted 1-dimensional representation The extracted 3D geometries are required to conduct studies on computational fluid dynamics, particle transfer and deposition, ventilation, stress analysis and deformation simulations. Marching cubes algorithm [73] is a very well established method implemented in FAST [32] allowing the generation of 3D geometric models from airway segmentation label maps. The constriction simulation method aims to generate 3D tubular surface structures with smaller diameters. To this end, Laplacian surface contraction offers a solution that deforms the geometry pushing the vertices towards the direction of the inward normals. The extension of the extracted centerline generates a predictive representation of the bronchial tree given the available space. However, for the outcome of space-filling algorithms to be useful in fluid dynamics simulation, particle deposition simulation or stress finite element analysis based studies, the boundary conditions in the form of triangular 3D meshes need to be defined. Initially, as a simplified approach, to define the diameter of each generation we can employ the power law consistent with Murray’s law of symmetric branching [33, 34]. (3) where d0 denotes the branch diameter of the trachea and dz the branch diameter for generation z. Furthermore, if we take into account that each branch demonstrates different branching angle and diameter properties, the relation between airway diameter (d) and branching angles (θ) is based on the following rules validated by Kamiya et al. [37] and Kitaoka et al. [2]: (4) (5) where the index 0 stands for the parent branch, and the indices 1 and 2 for the two children branches, respectively. To reconstruct the lung surface we employ a point cloud sampling approach as input for Poisson surface reconstruction. The outcome is a smooth surface with smooth transitions instead of abrupt transitions in the intersection with the original tubular meshes. The tubular-shaped point cloud is sampled using a uniform random distribution. A clean-up step, visualized in Fig 7, ensures that no point can be found in distance less than the prescribed diameter of every available branch. The resulting point cloud is used to compute normals. A bilateral normal smoothing [74] function prepares the point cloud for Poisson surface reconstruction [75]. smoothing the point normals. This step facilitates the surface reconstruction in bifurcations and transitional parts. Furthermore, since the directed graph is extracted where each point on the centerline corresponds to a point on the lung surface it is possible to further deform the surface with a custom function or pattern. The generated surface for seven and ten generations is presented in Fig 8. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 7. Outcome of point cloud generation process. Extra refinements remove the inner points facilitating normal estimation for Poisson surface reconstruction. https://doi.org/10.1371/journal.pone.0230259.g007 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 8. Generation of surface from the extended bronchial tree centerline. https://doi.org/10.1371/journal.pone.0230259.g008 2.5 Simulation of constrictive pulmonary diseases’ effect on airway tree This section aims to provide the methodology for simulation of broncho-constriction allowing to subsequently estimate the dynamic behaviour of the lung airways in the case of an exacerbation event. A bronchial tree 3D geometry is the input for this process yielding as output contracted airways. The proposed geometry contraction procedure is presented by Nousias et al. [45] and Lalas et al. [46] and is an extension of the work of Au et al. [76] facilitating a Laplacian smoothing process that shifts vertices along the estimated curvature normal direction. The airway geometric model consists of connected triangles forming the boundary conditions. Each triangular mesh can be described as where is the set of vertices, is the set of edges and is the set of faces constituting the 3D surface. Each vertex i can be represented as a point vi = (xi, yi, zi), ∀i = 1, 2, ⋯, N. For each face fi, ∀i = 1, 2, ⋯, l we denote the centroid (6) The outward unit normal to the face fi (located at the centroid mi) is calculated as : (7) where are the vertices corresponding to face fi. Given the curvature flow Laplacian operator, the product δ = LV approximates the inward curvature flow normals [77]. The motivation for employing the curvature flow Laplacian operator [78] on the mesh is that its output is not affected by mesh density. Specifically, (8) where Ai is the one-ring area, κi is the local curvature and ni is the inward curvature flow normal of the ith vertex. The positions of the vertices satisfying LV = 0 result in a zero volume string-like mesh and can be used to simulate mesh contraction. However, since such an optimization problem has more than one solutions, further constraints are required [76]. Introducing weighting matrices can smoothly move vertex positions towards the direction of the inward unit normal by iteratively solving the following minimization problem (9) where corresponds to the vertex positions before the contraction at each iteration. The weighting matrices WH and WL regulate the mesh contraction and mesh attraction, respectively. Initially, we set them to and WH = I, where is the identity matrix, k a double constant experimentally tuned to 10−3 and A the average face area of the model. Eq (9) can be expressed as (10) The analytical solution can be formulated as (11) where matrices A and b are given by (12) After each iteration t we update the contraction and inflation weights to be used in iteration t + 1 so that and , where and are the original and the current one-ring area respectively. The Laplacian matrix for iteration t + 1, Lt+1 is also recomputed. On these grounds, to simulate broncho-constriction we require to reduce the airway diameter to a predefined level of narrowing is reached. This level is defined by certain termination criteria [45, 46]. Thus, a metric is required that measures the diameter of the bronchi under process. To estimate the degree of contraction of the airway’s geometry after each iteration, we employ a shape diameter function (SDF) based scheme [79] implemented in [80] that evaluates the local volume based on the estimated local diameter assigned to each face of the mesh, also known as raw SDF values. Measuring the volume before and after the Laplacian contraction iteration can set the termination criteria. Fig 9 presents a simulation of constrictive pulmonary conditions. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 9. Simulation of constrictive pulmonary conditions. https://doi.org/10.1371/journal.pone.0230259.g009 3 Results 3.1 Dataset description For the evaluation of the aforementioned approaches we employed the dataset provided by VESSEL12 (VESsel SEgmentation in the Lung) challenge [81] and EXACT09 [17]. The VESSEL dataset is comprised of 20 anonymized scans in Meta (MHD/raw) format. The latter consists of 75 completely anonymized chest CT scans contributed by eight different institutions, acquired with several different CT scanner brands and models, using a variety of scanning protocols and reconstruction parameters. The conditions of the scanned subjects varied widely, ranging from healthy volunteers to patients showing severe abnormalities in the airways or lung parenchyma. Fig 6a to 6c present imaging instances of CT slices across the axial, coronal and sagittal planes. The generation of the initial airway surface, lung volume and 1-dimensional representation are performed using the FAST framework [32]. 3.2 Structural modelling and validation Our simulation framework processes the initial tree centerline and generates a structural estimation given the first three to four available generation and their morphometric characteristics i.e., lengths and diameters. To facilitate the comparison with morphometric data, we employed a publicly available dataset provided by Montesantos et al. [44] labelled as pone.0168026.s001. For the sake of self-completeness, the authors of [44] provided morphometric data extracted from HRCT images acquired at the University Hospital Southampton NHS Foundation Trust as a part of study described in [82, 83]. Data from seven healthy subjects and six patients with moderate or persistent asthma were included in the dataset. Asthmatic patients patients were diagnosed exacerbation-free for at least one month and were male non-smokers. A Sensation 64 slice HRCT scanner (Siemens, Enlargen, Germany) was utilized to capture 3D images from mouth to the base of the lungs. Subjects were posed in supine position and were instructed to perform slow exhalation. Groundtruth data for the development of bronchial tree models in [44] were extracted by Pulmonary Workstation 2 Software including 3 to 4 generations in the upper lobes and 6 to 7 generations in the lower lobes. For each subject, a morphology file includes the total lung volume of the subject lung (in cm3) and the percent volume per lobe while a translation file contains the airway connectivity, starting from the trachea to the terminal nodes. We used the generated trees from [44] to validate our approach and compare our results with relevant literature findings. Specifically, we compared the distributions of diameters, airway lengths and branching angles for each generation and the total number of airways for Horsfield and Strahler orders. In total 31204 acini were calculated being in agreement with the results reported by [6, 44]. Figs 10 and 11 present a comparison in terms of the number of airways for each level of Strahler and Horsfield orders. This comparison confirms that our model comes into agreement with pone.0168026.s001. Furthermore, distributions of airway lengths, branching angles and diameters were plotted for each generation, for AVATREE and pone.0168026.s001 [44]. Airway lengths maintain the same exponential decay pattern for both models. Differences appear in generations 1 to 4 that are distinctively defined by body size and anatomical features. The distribution of branching angles of our model is also confirmed by pone.0168026.s001 [44] maintaining a nearly linear decay with a very small rate. The distributions of diameters per generation are also observed to follow an exponential decay pattern. Both our model and pone.0168026.s001 [44] decay similarly after generation 4 validating the morphometric characteristics of the airway trees generated by our approach. Figs 12 to 14 present the distribution of airway length, branching angle and diameter for each generation for AVATREE and for pone.0168026.s001 [44]. Table 1 presents and overview of quantitative macroscopic figures for AVATREE and relevant studies. Branching ratios (RBH, RBS), diameter ratios (RDH, RDS) and length ratios RLH, RLS) were calculated for Strahler and Horsfield ratios denoted as *S and *H respectively. Specifically, RBH, RDH and RLH were calculated equal to RBH = 1.74, RDH = 1.259 and RLH = 1.26±1.01. Montesantos et al. [44] reported RBH = 1.56, RDH = 1.115 and RLH = 1.13 respectively. Additionally, RBS, RDS and RLS were calculated equal to RBS = 2.35, RDS = 1.25 and RLS = 1.23±1.02 and are close to the figures provides by relative studies [1, 44] as Table 1 reveals. Likewise, rate of decline for diameters per generation RD was calculated to RD = 0.83±0.21, being in agreement to [44]. Average branching angle θ for our model was calculated to 32.44±28.95 comparable to [44] reporting a θ equal to 42.1±21.4. Finally, Figs 15 and 16 present bronchial tree 1-dimensional representations extended up to 12 and 23 generations respectively. Additionally, for both generated models the surface has been reconstructed for the first 7 generations. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 1. Structural features comparison. https://doi.org/10.1371/journal.pone.0230259.t001 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 10. Number of airways for each Strahler order for our model, “AVATREE”, and “pone.0168026.s001”. https://doi.org/10.1371/journal.pone.0230259.g010 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 11. Number of airways for each Horsfield order for our model, “AVATREE”, and “pone.0168026.s001”. https://doi.org/10.1371/journal.pone.0230259.g011 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 12. Distribution of airway lengths for each generation for AVATREE and pone.0168026.s001. https://doi.org/10.1371/journal.pone.0230259.g012 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 13. Distribution of branching angles for each generation for AVATREE and pone.0168026.s001. https://doi.org/10.1371/journal.pone.0230259.g013 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 14. Distribution of diameters for each generation for AVATREE and pone.0168026.s001. https://doi.org/10.1371/journal.pone.0230259.g014 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 15. Estimation of bronchial tree for 12 generations. Surface reconstruction was performed for the first 7 generations. https://doi.org/10.1371/journal.pone.0230259.g015 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 16. Estimation of bronchial tree for 23 generations. Surface reconstruction was performed for the first 7 generations. https://doi.org/10.1371/journal.pone.0230259.g016 3.3 Open-source library & Graphical User Interface The presented components of AVATREE are provided as an open-source solution publicly available at (https://gitlab.com/LungModelling/avatree) accompanied by a graphical user interface (GUI). The implemented application programming interface (API) includes the following modules, a) input-output functionalities, b) 1-dimensional representation tools including centerline extraction, graph generation, derivation of graph node properties c) bronchial tree extension tools extending the 1-dimensional representation to the desired number of generations, d) 3D surface generation and processing tools and e) airway narrowing simulation tools. The GUI, presented in Fig 17, employs the set of functionalities defined by AVATREE and is comprised of four panels, namely data input and output panel, area selection panel, segmentation panel and broncho-constriction simulation panel. Through the GUI the user can load a 3D model, select the area to be processed, as Fig 18 visualizes, and use the narrowing functionalities to reduce the airway diameter by the desired percentage. The amount of narrowing depends on the number of iterations and contraction weight multiplier. In Fig 18 an airway of the first generation was constricted by 66%. The deformed surface introduced into computational fluid dynamics can provide insight into the breathing pattern and drug delivery in asthmatic lungs [84]. In the segmentation panel, the surface faces can be classified based on local properties. The one illustrates the shape diameter function (SDF) [79], while the other one the 3D surface according to the generation number. The results are visualized in Fig 19. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 17. User interface. https://doi.org/10.1371/journal.pone.0230259.g017 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 18. Broncho-constriction simulation. Airway of second generation narrowed at 34% of original diameter. https://doi.org/10.1371/journal.pone.0230259.g018 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 19. Surface annotation with a) SDF function visualizing local diameter b) airway generation. https://doi.org/10.1371/journal.pone.0230259.g019 3.1 Dataset description For the evaluation of the aforementioned approaches we employed the dataset provided by VESSEL12 (VESsel SEgmentation in the Lung) challenge [81] and EXACT09 [17]. The VESSEL dataset is comprised of 20 anonymized scans in Meta (MHD/raw) format. The latter consists of 75 completely anonymized chest CT scans contributed by eight different institutions, acquired with several different CT scanner brands and models, using a variety of scanning protocols and reconstruction parameters. The conditions of the scanned subjects varied widely, ranging from healthy volunteers to patients showing severe abnormalities in the airways or lung parenchyma. Fig 6a to 6c present imaging instances of CT slices across the axial, coronal and sagittal planes. The generation of the initial airway surface, lung volume and 1-dimensional representation are performed using the FAST framework [32]. 3.2 Structural modelling and validation Our simulation framework processes the initial tree centerline and generates a structural estimation given the first three to four available generation and their morphometric characteristics i.e., lengths and diameters. To facilitate the comparison with morphometric data, we employed a publicly available dataset provided by Montesantos et al. [44] labelled as pone.0168026.s001. For the sake of self-completeness, the authors of [44] provided morphometric data extracted from HRCT images acquired at the University Hospital Southampton NHS Foundation Trust as a part of study described in [82, 83]. Data from seven healthy subjects and six patients with moderate or persistent asthma were included in the dataset. Asthmatic patients patients were diagnosed exacerbation-free for at least one month and were male non-smokers. A Sensation 64 slice HRCT scanner (Siemens, Enlargen, Germany) was utilized to capture 3D images from mouth to the base of the lungs. Subjects were posed in supine position and were instructed to perform slow exhalation. Groundtruth data for the development of bronchial tree models in [44] were extracted by Pulmonary Workstation 2 Software including 3 to 4 generations in the upper lobes and 6 to 7 generations in the lower lobes. For each subject, a morphology file includes the total lung volume of the subject lung (in cm3) and the percent volume per lobe while a translation file contains the airway connectivity, starting from the trachea to the terminal nodes. We used the generated trees from [44] to validate our approach and compare our results with relevant literature findings. Specifically, we compared the distributions of diameters, airway lengths and branching angles for each generation and the total number of airways for Horsfield and Strahler orders. In total 31204 acini were calculated being in agreement with the results reported by [6, 44]. Figs 10 and 11 present a comparison in terms of the number of airways for each level of Strahler and Horsfield orders. This comparison confirms that our model comes into agreement with pone.0168026.s001. Furthermore, distributions of airway lengths, branching angles and diameters were plotted for each generation, for AVATREE and pone.0168026.s001 [44]. Airway lengths maintain the same exponential decay pattern for both models. Differences appear in generations 1 to 4 that are distinctively defined by body size and anatomical features. The distribution of branching angles of our model is also confirmed by pone.0168026.s001 [44] maintaining a nearly linear decay with a very small rate. The distributions of diameters per generation are also observed to follow an exponential decay pattern. Both our model and pone.0168026.s001 [44] decay similarly after generation 4 validating the morphometric characteristics of the airway trees generated by our approach. Figs 12 to 14 present the distribution of airway length, branching angle and diameter for each generation for AVATREE and for pone.0168026.s001 [44]. Table 1 presents and overview of quantitative macroscopic figures for AVATREE and relevant studies. Branching ratios (RBH, RBS), diameter ratios (RDH, RDS) and length ratios RLH, RLS) were calculated for Strahler and Horsfield ratios denoted as *S and *H respectively. Specifically, RBH, RDH and RLH were calculated equal to RBH = 1.74, RDH = 1.259 and RLH = 1.26±1.01. Montesantos et al. [44] reported RBH = 1.56, RDH = 1.115 and RLH = 1.13 respectively. Additionally, RBS, RDS and RLS were calculated equal to RBS = 2.35, RDS = 1.25 and RLS = 1.23±1.02 and are close to the figures provides by relative studies [1, 44] as Table 1 reveals. Likewise, rate of decline for diameters per generation RD was calculated to RD = 0.83±0.21, being in agreement to [44]. Average branching angle θ for our model was calculated to 32.44±28.95 comparable to [44] reporting a θ equal to 42.1±21.4. Finally, Figs 15 and 16 present bronchial tree 1-dimensional representations extended up to 12 and 23 generations respectively. Additionally, for both generated models the surface has been reconstructed for the first 7 generations. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 1. Structural features comparison. https://doi.org/10.1371/journal.pone.0230259.t001 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 10. Number of airways for each Strahler order for our model, “AVATREE”, and “pone.0168026.s001”. https://doi.org/10.1371/journal.pone.0230259.g010 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 11. Number of airways for each Horsfield order for our model, “AVATREE”, and “pone.0168026.s001”. https://doi.org/10.1371/journal.pone.0230259.g011 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 12. Distribution of airway lengths for each generation for AVATREE and pone.0168026.s001. https://doi.org/10.1371/journal.pone.0230259.g012 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 13. Distribution of branching angles for each generation for AVATREE and pone.0168026.s001. https://doi.org/10.1371/journal.pone.0230259.g013 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 14. Distribution of diameters for each generation for AVATREE and pone.0168026.s001. https://doi.org/10.1371/journal.pone.0230259.g014 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 15. Estimation of bronchial tree for 12 generations. Surface reconstruction was performed for the first 7 generations. https://doi.org/10.1371/journal.pone.0230259.g015 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 16. Estimation of bronchial tree for 23 generations. Surface reconstruction was performed for the first 7 generations. https://doi.org/10.1371/journal.pone.0230259.g016 3.3 Open-source library & Graphical User Interface The presented components of AVATREE are provided as an open-source solution publicly available at (https://gitlab.com/LungModelling/avatree) accompanied by a graphical user interface (GUI). The implemented application programming interface (API) includes the following modules, a) input-output functionalities, b) 1-dimensional representation tools including centerline extraction, graph generation, derivation of graph node properties c) bronchial tree extension tools extending the 1-dimensional representation to the desired number of generations, d) 3D surface generation and processing tools and e) airway narrowing simulation tools. The GUI, presented in Fig 17, employs the set of functionalities defined by AVATREE and is comprised of four panels, namely data input and output panel, area selection panel, segmentation panel and broncho-constriction simulation panel. Through the GUI the user can load a 3D model, select the area to be processed, as Fig 18 visualizes, and use the narrowing functionalities to reduce the airway diameter by the desired percentage. The amount of narrowing depends on the number of iterations and contraction weight multiplier. In Fig 18 an airway of the first generation was constricted by 66%. The deformed surface introduced into computational fluid dynamics can provide insight into the breathing pattern and drug delivery in asthmatic lungs [84]. In the segmentation panel, the surface faces can be classified based on local properties. The one illustrates the shape diameter function (SDF) [79], while the other one the 3D surface according to the generation number. The results are visualized in Fig 19. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 17. User interface. https://doi.org/10.1371/journal.pone.0230259.g017 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 18. Broncho-constriction simulation. Airway of second generation narrowed at 34% of original diameter. https://doi.org/10.1371/journal.pone.0230259.g018 Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 19. Surface annotation with a) SDF function visualizing local diameter b) airway generation. https://doi.org/10.1371/journal.pone.0230259.g019 4 Conclusion In this paper, we presented AVATREE, an end-to-end approach modeling the subject-specific airway tree that defines the personalized boundary conditions required for the simulation of pulmonary function. This particular personalization aspect refers to the extraction of the main airways and the lung volumes from imaging allowing the simulation of a personalized extended bronchial tree. The utter goal in this category of studies is to eventually predict gas flow and particle distribution in healthy and constricted bronchial trees. Modelling lung ventilation patterns can provide grounds for performance and fatigue estimation in high-frequency ventilation cases and give insight into drug delivery or even transfer of harmful particulates. Specifically, this work presents an open-source simulation framework that utilizes imaging data to provide patient-specific representations of the structural models of the bronchial tree. The extended 1-dimensional centerline facilitates the generic estimation of pulmonary function through lumped 0-dimensional studies and allows the generation of probabilistic confidence maps of airway generation data. Such a visualization could be exploited by airway tree segmentation methodologies to improve the results further constraining the 3D space to be searched. Further elaborating on the benefits of the presented methodology, the generation of extended bronchial tree surface allows the assessment of airway functionality. Several studies available in the literature have employed computational fluid dynamics to predict flow and particle transfer patterns inside the conducting regions of the bronchial tree. Generating the 3D mesh constituting the surface defines the boundary conditions for this category of studies. Furthermore, surface deformation functionalities allow simulating broncho-constriction, which is the main feature in constrictive pulmonary diseases such as asthma and COPD. Existing approaches of medication adherence in asthma and COPD patients are usually based on the analysis of breathing signals acquired with acoustic sensors attached on inhaler devices [85] [86]. The concept of such studies is to facilitate self-management by guiding the patients to improve their inhaler usage technique [87]. AVATREE could contribute to this type of analysis by estimating the effectiveness of inhaled medication based on personalized imaging data and particle deposition simulations [46]. Both automated and UI guided solutions are provided by the presented open-source solution enabling users to simulate pathological conditions in asthmatic patients guided by imaging priors data from healthy subjects.
TI - AVATREE: An open-source computational modelling framework modelling Anatomically Valid Airway TREE conformations
JF - PLoS ONE
DO - 10.1371/journal.pone.0230259
DA - 2020-04-03
UR - https://www.deepdyve.com/lp/public-library-of-science-plos-journal/avatree-an-open-source-computational-modelling-framework-modelling-oxhZPdDbjW
SP - e0230259
VL - 15
IS - 4
DP - DeepDyve
ER -