Computers in Biology and Medicine 62 (2015) 25–32
Contents lists available at ScienceDirect
Computers in Biology and Medicine journal homepage: www.elsevier.com/locate/cbm
Mathematical model of a heterogeneous pulmonary acinus structure Kenichiro Koshiyama n, Shigeo Wada Graduate School of Engineering Science, Osaka University, Japan
art ic l e i nf o
a b s t r a c t
Article history: Received 30 October 2014 Accepted 31 March 2015
The pulmonary acinus is a gas exchange unit distal to the terminal bronchioles. A model of its structure is important for the computational investigation of mechanical phenomena at the acinus level. We propose a mathematical model of a heterogeneous acinus structure composed of alveoli of irregular sizes, shapes, and locations. The alveoli coalesce into an intricately branched ductal tree, which meets the space-filling requirement of the acinus structure. Our model uses Voronoi tessellation to generate an assemblage of the alveolar or ductal airspace, and Delaunay tessellation and simulated annealing for the ductal tree structure. The modeling condition is based on average acinar and alveolar volume characteristics from published experimental information. By applying this modeling technique to the acinus of healthy mature rats, we demonstrate that the proposed acinus structure model reproduces the available experimental information. In the model, the shape and size of alveoli and the length, generation, tortuosity, and branching angle of the ductal paths are distributed in several ranges. This approach provides a platform for investigating the heterogeneous nature of the acinus structure and its relationship with mechanical phenomena at the acinus level. & 2015 Elsevier Ltd. All rights reserved.
Keywords: Pulmonary biomechanics Biological microstructure Pulmonary physiology Gas exchange Computational geometry Structure modeling Rat acinus
1. Introduction The pulmonary acinus is a gas exchange unit that is located distal to the terminal bronchioles and consists of a complex microstructure of ducts and alveoli [1–3]. The mechanical phenomena at the acinus level, i.e., acinar flow, gas diffusion, particle depositions, or mechanical responses of lung parenchyma to stretches, underlie the physiological behavior of the lungs. A detailed understanding of the pulmonary acinus is required for recent developments in medical engineering techniques, e.g., pulmonary drug delivery [4,5] or mechanical ventilation protocols [6,7]. To date, pulmonary researchers have investigated acinus phenomena using theoretical and numerical approaches [8–15]. The structure model of the acinus (i.e., the shape and size of the alveoli and the patterns of the ductal tree) affects and restricts the numerical results, and hence the choice of structure model is important. Various image-based [16,17] and mathematical models [11,18–21] have been proposed. Image-based models represent part of the real acinus structure, but are sensitive to the image resolution and segmentation algorithm. Hence, their application to numerical analysis is still restricted. In contrast, mathematical models can easily mimic general characteristics of the acinus
n Correspondence to: Graduate School of Engineering Science, Osaka University, Toyonaka, 560-8531, Japan. Tel.: þ 81 6 6850 6174; fax: þ81 6 6850 6173. E-mail address: [email protected] (K. Koshiyama).
http://dx.doi.org/10.1016/j.compbiomed.2015.03.032 0010-4825/& 2015 Elsevier Ltd. All rights reserved.
structure. Two types of models are frequently used for further numerical analysis. The first type is composed of straight tubes that enable regular sphere caps or polyhedra to coalesce on the side wall (i.e., duct-based model) and is typically used to analyze acinar flows [21]. The second type of model is formed from polyhedra of regular size and shape (i.e., polyhedra-based model) and is used to analyze the mechanical response of lung parenchyma to stretching [22]. Conventional models rely on the average morphological characteristics of the acinus structure. Consequently, they are too basic to capture the heterogeneous nature of the alveolar and ductal patterns observed in real acinus structures. Computations using simplified acinus models provide approximated features of mechanical phenomena, but allow only a limited investigation of phenomena related to the heterogeneous nature of the acinus structure. In fact, variations of mechanical phenomena in the lungs (i.e., regional differences [23]) have long been recognized. Advances in mechanobiology mean that the detailed resolution of the strain field in alveolar walls is of great interest to researchers in lung cell mechanics [24]. Considering lung parenchyma mechanics, the pathway for stress transmission from organ to cells or cytoskeletons and vice versa is through the extracellular matrix. The parenchyma architecture i.e., detailed acinus structure, is one of important factors affecting such mechanotransduction [13,25]. The heterogeneity of the acinus structure on acinar flow remains important in biofluidics [21]. Furthermore, the acinus structure segmented from high-resolution CT images has been analyzed to micrometer resolution [26–28], which has the potential to provide
26
K. Koshiyama, S. Wada / Computers in Biology and Medicine 62 (2015) 25–32
quantitative data about the heterogeneity of the acinus structure. Thus, to understand the heterogeneous nature of the acinus structure and its relationship to mechanical phenomena in the lungs, it is vital to develop a parametric heterogeneous acinus structure model and investigate the potential and limitations of the modeling technique. In this paper, we propose a method for generating a heterogeneous acinus structure model composed of alveoli of irregular sizes, shapes, and locations that coalesce into an intricately branched ductal tree. Instead of including all the details of the actual acinus structure, we focus on the gas exchange unit of the 3D space-filling structure, represented by an assemblage of polyhedra. To our knowledge, several groups have developed heterogeneous acinus structure models [29,30], but they have focused on the heterogeneity of the ductal tree. A distinctive feature of our model is that we parametrically and explicitly include the distributions of both alveolar and ductal structures, without significantly departing from the typical polyhedra-based acinus model [18,19]. This is achieved by introducing Voronoi and Delaunay tessellations [31] and simulated annealing optimization [32,33]. The abundance of quantitative and qualitative morphological data about the acinus structure [2,34] allow us to verify our approach by modeling the acinus structure of a rat and comparing the model output with previous results. Our paper is organized as follows. In Section 2, we describe the basic methodology of heterogeneous modeling, present the model setting for the rat acinus structure, and introduce the structural parameters analyzed. Section 3 describes the sensitivity of the average structural parameters to the modeling method, distributions and variations of alveolar and ductal structures, segmentation of the 3D structure, model predictions, and limitations. Finally, Section 4 summarizes the paper.
2. Methods The approach presented here is a unified version of the polyhedra-based acinus structure modeling first presented by Denny and Schoroter [18]. The polyhedra-based modeling regards each polyhedron as an airspace region and the faces as tissue regions, i.e., alveolar or ductal walls. The approach is composed of two modeling steps: the first generates an assemblage of polyhedra, and the second generates a ductal tree by eliminating the common walls between these polyhedra based on the optimized ductal tree when the gas exchange efficiency in the acinus is considered as the cost for the optimization. In the following, we describe the extension of this approach to express the heterogeneity of the acinus structure. We illustrate the 2D case for simplicity (Fig. 1), although the method can be extended to three dimensions.
that when the seed coordinates are set to the body-centered cubic lattice in a 3D space, the structure of the resultant Voronoi polyhedron becomes a truncated octahedron, which is a typical shape for alveolar or ductal airspaces [18,19,29]. 2.2. An intricately branched ductal tree A tree structure is generated by connecting the neighboring seeds. To consistently obtain the neighboring seeds in the heterogeneous system, we introduce the Delaunay tessellation (Fig. 1C) from the mathematical background of the Voronoi tessellation [31]. To generate a fully isolated structure, we do not consider connections to seeds in polygons whose edges open outward (boundary polygons) (Fig. 1B). A seed near the center of the top edge of a given region is designated as the entrance seed, and subsequent connections are randomly selected from the neighboring seeds (Fig. 1D). The tree generation obeys two geometrical constraints based on experimental findings from an acinar ductal tree [2,18]: (i) each seed has a unique path to the entrance seed (i. e., no cycle); and (ii) the number of daughter branches from a seed is less than 3, whereas the number of branches to the terminal seed, i.e., to the alveoli, is unlimited. Because the acinus is a gas exchange unit, we control the geometrical components related to gas exchange efficiency by optimizing the tree structure based on the following criteria. The seeds are numbered, and a solution in the optimization provides each seed with information about one parent (upstream) seed number and the daughter (downstream) seed numbers (i.e., connective information). The neighbor of a solution is obtained by randomly changing the parent seed of a randomlyselected seed. As the gas exchange in real acini involves a complex set of factors (e.g., gas diffusion and convection in airways, gas exchange between airways and surrounding lung capillary networks, blood circulation or respiration conditions [35]), the selection of cost for the optimization is arbitrary. We assume that efficient gas exchange in acini is related to a shortening of ductal path lengths for gas diffusion and an increase in the number of alveoli for the region of the gas exchange which can be calculated from connective information with less computational costs. The same criterion has been used previously [18,34]. Accordingly, the cost to be minimized during optimization includes the cost C 1 of maximizing the number of alveoli (i.e., to minimize the number of ducts) C 1 ¼ Nt Na ;
where N t is the total number of seeds and N a is the total number of terminal seeds, as calculated from the connective information; and the cost C 2 of minimizing the average path length from the entrance to each terminal seed
2.1. An assemblage of airspace regions of irregular shapes and sizes C2 ¼ Firstly, seeds of polygons (polyhedra in 3D) are scattered in a given (square) space (Fig. 1(A)). The spatial coordinates are given by sets of pseudo-random real numbers taken from a uniform distribution between 0 to the distance in one direction of the space. We set the number density of seeds and the minimum distance as the conditions for seed generation. Computationally, for a given seed density, decreasing the minimum distance increases free space to generate seed coordinates, hence allowing various seed distributions. Polygons filling a given space are generated based on seeds by the digitized Voronoi tessellation technique, which divides a region by allocating points to the closest seeds [31] (Fig. 1A and B). Accordingly, each seed has information about edges of its own polygon (faces of its own polyhedron in 3D). It should be noted
ð1Þ
Na 1 X l; Na i ¼ 1 i
ð2Þ
where li is the path length to the ith terminal seed, which is the sum of the lengths between the connected seeds' coordinates in the system. The total cost C is defined as C ¼ ω C 1' þ ð1 ωÞC 2' ;
ð3Þ
where ω is the weight of each cost component, 0 r ω r 1, and the apostrophe denotes that each cost is normalized according to its values at ω ¼ 0 or 1. For the optimization, we introduce a simulated annealing algorithm, where the optimized tree structure is searched by modifying the connections between seeds and comparing the costs before and after the modification. Specifically, the parent of
K. Koshiyama, S. Wada / Computers in Biology and Medicine 62 (2015) 25–32
27
Fig. 1. Schematic diagram of heterogeneous acinus structure modeling: random seed distribution (A), Voronoi diagram (B), Delaunay diagram (C), initial ductal tree structure (D), optimized ductal tree structure (E), and resulting acinus structure model (airspace is shown in gray) (F).
Table 1 Rat acinus modeling parameters and their reproducibility. Dmin [μm]
Lx [mm]
ρ Vt F [-/mm3] [106 μm3] [%]
Va [μm3]
RV t [%] RF [%] RVa [%]
50 40 30 Experiment [2]
0.97 0.98 0.99 –
5478 5312 5153 –
196 190 187 194
99 99 98 100
527 527 532 522
76 76 75 80
95 94 94 100
99 98 98 100
Dmin : minimum distance between seeds; Lx : distance in one direction of the cubic system; ρ: number density of seeds; V t : average total airspace volume; F: average fraction of airspace; and V a : average alveolar volume. RVt , RF , and RV a denote the reproducibility of V t , F, and V a (see Section 2.4).
a randomly chosen seed is changed to a neighboring seed on the Delaunay sides, and then the new cost is calculated via Eq. (3) with the new tree structure. When the new cost is smaller than the previous one, the new tree structure is accepted; even when the new cost is larger than the previous one, it may be accepted on statistical terms (i.e., under the Metropolis criterion). The probability of acceptance of a larger cost decreases with the number of modifications, and the modification is iterated until all costs have converged. After the optimized tree structure was obtained (Fig. 1E), the polygonal edges pierced by the optimized tree are eliminated. As each seed has information about edges of its own polygon, eliminating the common edges of connecting seeds results in all the polygons in the system being vented from the entrance (Fig. 1F). 2.3. 3D acinus model settings To visualize and analyze the 3D geometry of the acinus model, a volume rendering technique is required. In practice, this is achieved by forming the 3D digitized Voronoi polyhedra using
an object of arbitrary shape composed of a number of fine cubic building blocks (voxels). The dimensions of the object, voxel size, number of voxels representing the acinus tissue (i.e., alveolar wall thickness), and number of seeds are related to the acinar and alveolar sizes. These parameters are set so that the total airspace volume, average alveolar volume, and volume density of the airspace in the model match those determined experimentally. The three volume-related data are the modeling condition for 3D acinus modeling. For simplicity, we assume that the acinus is cubic, all alveoli belong to a single acinus, and the number of seeds, the voxel size, and the wall thickness are constant. Furthermore, the seed coordinates are given by a pseudo-random number generator with some minimum distance between seeds, by which we attempt to control spatial distributions of seeds. For the rat acinus structure, the number of seeds, the voxel size, and wall thickness are set to 5000, 2 μm, and 4 μm, respectively; other modeling parameters are listed in Table 1. For generating irregular seed distribution, sets of pseudo-random numbers, taken from a uniform distribution in the range ð0; 1Þ, are scaled by the distance in one direction of the cubic system and are used as spatial coordinates of seeds. The minimum distances are 30, 40, and 50 μm. Note that, in practice, the modeling parameters for each minimum distance are set empirically to fulfill the volumerelated data from experiments. Accordingly, the resultant number densities of seeds for different minimum distances differ slightly (Table 1), but controlling seed distributions with minimum distance work with the given conditions and assumptions (see Section 3.2). #Although other structural information is available for the rat acinus (e.g., number of alveoli per acinus, average path length, etc. [2]), volume-related information is relatively reliable, as it can be obtained by simply counting the voxels in acinus images, whereas detailed alveolar or ductal structural information requires interpretation (e.g., skeletonization for the average path length, and estimation from small airway samples or 2D images for the
28
K. Koshiyama, S. Wada / Computers in Biology and Medicine 62 (2015) 25–32
structures increase with the weight, as expected, and exhibit a plateau region from 0:2 to 0:8, which corresponds to the previous regular polyhedra-based model (Fig. 4 in [34]). Additionally, the trend is insensitive to the initial seed distribution and the minimum distance between seeds (Fig. S2 in Supplementary information file A). The maximum number of generations, which has been shown to be important in verifying the path structure [34], and other path structures also exhibit a plateau region when the weight is not close to 0 or 1 (Fig. 2B and Fig. S3 in Supplementary information file A). These results suggest that, as long as we consider the two costs in Eq. (3) for gas exchange efficiency, the average path structure is essentially invariant. This conclusion is confirmed by the statistics of experimental results [2] and the previous regular structure model [34]. Although this insensitivity to weight (away from the extreme values) is evident in the current and previous [34] polyhedrabased models, there are slight differences in the invariant values. The current model gives approximately 2800 alveoli in a maximum of 14 generations, whereas the previous model suggests around 2200 alveoli in 12 generations. This difference can be explained by the explicit modeling of the wall thickness and alveolar size distributions in the current model, as the total number of polyhedra ( 3200 after eliminating boundary seeds) for the volume-related modeling condition (Table 1) is larger than that of the previous model ( 2700) [34].
number [27]). Therefore, in the current study, we use volumerelated data as the modeling conditions. For the simulated annealing, we introduce the Johnson–Aragon–McGeoch–Schevon algorithm [32] (detailed algorithm and parameters are given in Supplementary information file A). To verify the effects of the weight in Eq. (3) on the resulting acinus structure, we change the weight from 0 to 1. For each weight, we analyze the data from three different initial seed distributions. The modeling and analysis use in-house code written in C and Fortran 90 with the aid of the Qhull computational geometry algorithms [36] and the Mersenne Twister random seed generator. The source code is available from the authors upon request. To visualize the 3D geometry, we render sets of sliced images and reconstruct them using the general visualization tool Amillas. 2.4. Analysis Based on the definition given in previous studies [1,34], we regard the path from a branching seed to the next branching seed as one generation, ignoring the path to the terminal seed (alveoli). The complexity of the ductal tree is characterized by the distributions of the tortuosity and of the branching angle. The tortuosity is estimated as the ratio between the actual path length and the linear distance from the entrance seed to the terminal branching seed, where the actual path length is the sum of the path lengths per generation. The branching angle is that between the direction to the periphery of the parent's path and that of the daughter's path. The reproducibility R of the experimental information is calculated as R ¼ 1 jV V vj, where v is the average value given by the model and V is the experimental value. R is a distance measure, i.e., absolute deviation, between the model value and the experimental value, normalized by the experimental value. In order to mean that, at R ¼ 1:0(100%), the model value is equivalent to the experimental value, the distance is shifted one unit to the positive. R is only valid for 0 r v r 2Vas R becomes negative for v 4 2V.
3.2. Structural distributions for ω ¼ 0:5 and effects of minimum distance Because of the insensitivity of the average structure to the cost weighting, it is worth investigating the distributions of various structural parameters for a particular value (ω ¼ 0:5). One important characteristic of the polyhedra-based acinus model is that it expresses certain path structure distributions, i.e., path length, number of generations, tortuosity, and branching angles [18] (see Fig. 3C–F). In the current model, we extend this feature to express the distributions of alveolar shapes and sizes (Fig. 3A and B). By changing the initial seed distribution, we can create structure models in which the shape of each distribution is different while the average values remain almost constant. Fig. 4 shows the effects of the minimum distance Dmin on the alveolar shape and volume. For all Dmin , the peak value of the number of polyhedra faces is 14, which corresponds to the number of faces in the truncated octahedron often used to give a mathematical representation of the airspace in acinus [19]. The distributions of the number of faces and volume of alveoli become narrower as the minimum distance increases. Interestingly, the path structure distributions are insensitive to changes in this
3. Results and discussion 3.1. Effect of weight on the acinus structure To verify the sensitivity of the optimized structure to the weighting parameter in the heterogeneous acinus modeling, we first examined the sensitivity of the number of alveoli, average path length, and maximum number of generations. Fig. 2A shows the effect of ω on the number of alveoli and the average path length, which are directly optimized via costs C 1 and C 2 . Both
2500
3 No. of alveoli Av. path length
2000
2
1500
1
1000
0
0.2
0.4
0.6
ω[−]
0.8
1
0
Number of generations [ - ]
4
Average path length [ mm ]
Number of alveoli [ - ]
20 3000
18 16 14 12 10 0
0.2
0.4
0.6
0.8
1
ω[−]
Fig. 2. Effect of cost weight on the number of alveoli and average path length (A) and the maximum number of generations (B) for Dmin ¼ 40. Average values from three different seed distributions are shown. Error bars in Fig. 2B denote the standard deviation.
K. Koshiyama, S. Wada / Computers in Biology and Medicine 62 (2015) 25–32
30
30 20
30 20
Model 2
10
Model 1
10
Frequency [%]
Frequency [%]
20
Model 1
10
30 20
30 20
Model 2
10 30 20
Model 3
10 0
0
5
10
15
20
25
Model 3
10 0
30
0
3.0
4.0 2 ×10
30
20
20
Model 1
Model 1
10
Frequency [%]
10
Frequency [%]
2.0
Alveolar volume [ µm ]
30
30 20
Model 2
10 30 20
30 20
Model 2
10 30 20
Model 3
10
Model 3
10 0
0.5
1.0
1.5
0
2.0
0
Path length [ mm ]
10
15
20
30
20
20
Model 1
30 20
Model 2
10
Model 1
10
Frequency [%]
10
30
30 20
Model 2
10 30 20
20
Model 3
10 0 1.0
5
Number of generations [ - ]
30
Frequency [%]
1.0
3
Number of faces [ - ]
0
29
1.5
2.0
2.5
Model 3
10 3.0
Tortuosity [ - ]
0
0
45
90
135
180
Angle [ - ]
Fig. 3. Representative results about distributions of the number of faces (A), alveolar volume (B), path length (C), number of generations (D), tortuosity (E), and path angle (F) for Dmin ¼ 40 with ω ¼ 0:5. Model results generated from three different initial seed distributions are labeled as Model 1, 2, and 3.
minimum distance (Fig. S4 in Supplementary information file A). Because we have used a cubic acinus model with almost the same total volume for different Dmin (Table 1), we suspect that the shape of the acinus might affect the path structure distributions. We are now developing an ellipsoidal acinus model, and are analyzing the shape and volume information of individual acini from CT images [37]. The effects of acinar shape on the path structure will be investigated in more detail in the near future. 3.3. Variations of acinus structure during optimization The acinus models presented in Sections 3.1 and 3.2 and in previous studies [18,34] are based on the optimized tree structure. During optimization, the tree experiences various structures.
Fig. 5A presents representative results for the variation in average path length and number of alveoli. The two structural parameters change on each iteration to fulfill the modeling conditions (Table 1). Depending on the weight and initial seed distribution, the paths of the variations change (Fig. S5). For the modeling conditions of volume-related experimental data (Table 1), the average path length and the number of nonterminal seeds are minimized regardless of experimental values. As experimental values for the number of alveoli and average path length are available ( 1542 and 1.287 mm, respectively [2]), the reproducibility of these parameters can be examined along the path (Fig. 5B). Interestingly, the reproducibility approaches 1, before moving away during optimization. This indicates that, in the given condition, the actual rat acinus structure is not fully
30
K. Koshiyama, S. Wada / Computers in Biology and Medicine 62 (2015) 25–32
30
40 30
20 Dmin = 30 mm
Frequency [ % ]
Frequency [ % ]
10 30 20 Dmin = 40 mm 10 30 20
0
0
5
15
10
20
25
30 Dmin = 40 mm
20 10 40 30 20 10
Dmin = 50 mm
10
Dmin = 30 mm
20 10 40
30
0
Dmin = 50 mm
0
1.0
2.0
3.0
4.0
Alveolar volume [ µm 3 ]
Number of faces [ - ]
5.0
×10 2
Reproducibility
Average path length [ mm ]
Fig. 4. Effects of the minimum distance between seeds on the distributions of the number of faces (A) and alveolar volume (B). Distributions were calculated using data obtained from three models for each Dmin under ω ¼ 0:5.
6 5 4 3 2 1 0 1000
the distributions of alveolar and ductal structural parameters (Figs. 2–5) mean that the current model provides a platform for computational investigations of the heterogeneous nature of the acinus structure and its relationships with mechanical phenomena at the acinus level.
Initial tree Optimized tree
1500
2000 Number of alveoli [ - ]
2500
3000
1.0 Average path length
0.8 0.6 0.4 0.2 0
Number of alveoli Total
10 5 10 6 Number of iterations [ - ]
107
Fig. 5. Representative relationship between average path length and number of alveoli during optimization for ω ¼ 0:6 (A), and the change in reproducibility (B). In (B), the black solid line represents the reproducibility of the number of alveoli, gray dashed line represents average path length, and black bold line denotes total reproducibility.
optimized for gas exchange efficiency, at least with the given cost functions and the modeling assumptions in Section 2, and other mechanical factors may be involved in the morphology.
3.4. Segmentation of heterogeneous acinus structure models Segmentations of the airspace or tissue region are valuable for further numerical calculations. For this purpose, we explicitly modeled the alveolar wall and performed an image rendering (Fig. 6). After determining the ductal tree, and eliminating the walls that the ductal tree pierces, in a 3D cube (Fig. 6A and B), slices of the binarized images (in which tissue voxels and airspace voxels are shown in black and white, respectively) are stacked (Fig. 6C). This is a similar approach to the generation of 3D models from 2D images [2,28]. When the black voxels are reconstructed by general image processing tools (Fig. 6D), the model can be readily rendered as a 3D mesh for mechanical analyses, e.g., the mechanical response of the parenchyma to stretches, and when the white voxels are reconstructed and meshed (Fig. 6E), the model can be used to analyze gas diffusion or acinar flow. Part of the acinus (alveolar ducts or alveolar sac) can also be reconstructed (Fig. 6F). The variety of the average structure (Fig. S5) and
3.5. Model predictions and limitations A distinctive feature of our model is that we parametrically and explicitly include the distributions of both alveolar and ductal structures (Figs. 3–6). In addition to the comparison with the average data (Table 1 and Fig. 5B), we compare the distributions (Figs. 3 and 4) with available experimental data. Hansen and Ampaya [38] analyzed airspace regions of human acini and presented six shapes for alveolar airspace, assuming symmetrical shapes with curved surfaces. Fung [19] proposed 14hedron is a simple alveolar shape. Although we did not assume symmetrical shapes, the peak value of the number of polyhedral faces was 14 and they were distributed from 8 to 25, depending on the minimum distance (Fig. 4A). Mercer and Crapo [2] reported that the alveolar volume is 194 721 μm3 (average 7 standard error) for similar acinar conditions used for our modeling (Table 1). Although they provided only the standard error, the standard deviation might be large, considering the number of alveoli they analyzed (sixty alveoli from different lungs). The distributions of alveolar volume and shape can be controlled in our modeling by changing the minimum distance between seeds (Fig. 4), while a regular model provides only airspace region of regular volume and shape [18]. Among the pathway structure distributions, we verify the length distribution (Fig. 3C) because the frequency distributions for four individual acini have been reported by Mercer and Crapo [2]. From their experiment, regardless of acini analyzed, the lengths range from 0.2 to 2.5 mm and the maximum relative frequency is under 30%. The similar range of the path length for average acinus has been presented by Rodriguez, et al. [3]. Likewise, the lengths from our model are within the range of the above experimental data (Fig. 3C). In addition, the experimental data about different acini suggest that the shape of the distribution depends on individual acini [2]. The difference of the distribution shape can be obtained in our models for different seed distributions (Fig. 3C). The variety of alveolar and pathway structures in acinus has been widely observed in experiments [1–3,27,38–42]. The average data and their distributions depend on, species [41], inflation level [27,40], region in lungs [23], pathological conditions [42], or ages
K. Koshiyama, S. Wada / Computers in Biology and Medicine 62 (2015) 25–32
31
Fig. 6. Representative 3D heterogeneous acinus rendering: seeds distributed in a cube (A), optimized path tree (B), stack of 2D images (C), wall region rendering (D), airspace rendering (E), and part of the acinus model (F). See the text for more details.
[28]. Strictly, methodological issues (fixation of samples, imaging devices, segmentation algorithm, etc.) may also affect the data because of the delicate microstructure of acinus. Hence, great care is needed for the quantitative comparison of model predictions with experimental data. Nevertheless, as shown in Figs. 3–6, the architecture of acinus model presented here includes distributions of, not only pathway structures, but also alveolar shape, size, and locations, while the regular model has the distribution of pathway structures and alveolar location only [18]. Furthermore, we can parametrically control the distributions alveolar shape and size by setting minimum distance between seeds (Fig. 4). Besides the advantage, our approach has some limitations. First, the duct size may vary with the generation of ducts in the actual acinus [2,3], whereas our model assumes the size is independent of generation. Second, the alveolar shape may be restricted in the actual acinus [38], partly due to mechanical balance in acinus, whereas it is unlimited in our model. The first limitation is common in polyhedra-based models, and restricts our modeling to acini that do not include apparent respiratory bronchioles. Applied to human acinus, our model may correspond to the 1/8 subacinus [1], and hence, volume-related data should be taken from this region. The second point may affect the reproducibility (Fig. 5B). Because our model is a static model and does not involve any mechanics, further improvements should be made to e.g., the cost functions used in the current model. Alternatively, it may be interesting to investigate the effects of mechanics of alveolar walls (stiffness to the septal walls, negative pressure boundary, surface tension effects, etc. [25]) on the shape and size distributions by mechanical analysis using our structure models [8,16,25]. The assumptions used in the modeling, e.g., cubic acinar shape, constant wall thickness and number of seeds for setting parameters, etc., are also the limitations and should be addressed in the future studies. The resolution of in situ CT acinus imaging has reached the scale of the alveolar wall in various inflation levels [17,27], and our group has developed a semiautomatic segmentation algorithm that has the potential to provide the morphology of individual acini and a cluster of acini stemming from the same conducting airway [37]. Even in such detailed images, identifying ducts or alveoli is still challenging, but more sophisticated statistical data about alveolar and ductal distributions of individual acini will facilitate and complement our modeling, and vice versa.
4. Summary We have demonstrated that the combined use of Voronoi and Delaunay tessellations and the optimization of the ductal tree with simulated annealing can parametrically model a heterogeneous acinus structure composed of alveoli of irregular sizes, shapes, and locations that coalesce into an intricately branched ductal tree. The modeling condition considered average acinar and alveolar volume characteristics from published experimental information on rat acinus. In the model, the shape and size of alveoli and the length, generation, tortuosity, and branching angle of the ductal paths were distributed in several ranges. The distribution of the shape and size of the alveoli depends on the minimum distance between airspace regions in acinus models, while those of the path structure were found to be insensitive in cubic acinus models. The optimized tree structure did not fully reproduce the experimental values of the number of alveoli and average path length, indicating that the actual rat acinus structure is not fully optimized for gas exchange efficiency and that other mechanical factor may be involved in the morphology. The proposed technique can model the tissue and airspace of acinus without significantly departing from the typical acinus model used for various mechanical analyses. Therefore, the heterogeneous acinus model developed here can be used to investigate the heterogeneous nature of the acinus structure and its relationships with mechanical phenomena at the acinus level using computational approaches.
Conflict of interest statement The authors have no conflict of interests to declare with regard to this study.
Acknowledgments We would especially like to thank Prof. R.C. Schroter and Dr. T. Sera for their useful comments and suggestions. We also thank M. Yasui, A. Mizuno, T. Nishimura, and L. Xiao for their technical assistance. This work was partly supported by the Japanese Ministry of Education, Culture, Sports, Science and Technology (24240073).
32
K. Koshiyama, S. Wada / Computers in Biology and Medicine 62 (2015) 25–32
Appendix A. Supporting information Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j.compbiomed.2015. 03.032. References [1] B. Haefeli-Bleuer, E.R. Weibel, Morphometry of the human pulmonary acinus, Anat. Rec. 220 (4) (1988) 401–414. [2] R.R. Mercer, J.D. Crapo, Three-dimensional reconstruction of the rat acinus, J. Appl. Physiol. 63 (2) (1987) 785–794. [3] M. Rodriguez, et al., Pulmonary acinus—geometry and morphometry of the peripheral airway system in rat and rabbit, Am. J. Anat. 180 (2) (1987) 143–155. [4] N.R. Labiris, M.B. Dolovich, Pulmonary drug delivery. Part I: physiological factors affecting therapeutic effectiveness of aerosolized medications, Br. J. Clin. Pharmacol. 56 (6) (2003) 588–599. [5] J.S. Patton, J. Bukar, S. Nagarajan, Inhaled insulin, Adv. Drug Deliv. Rev. 35 (2–3) (1999) 235–247. [6] K.W. Delank, et al., Diagnosis and treatment modalities of sinunasal inverted papillomas, Laryngo-Rhino-Otol. 79 (4) (2000) 226–232. [7] J. Pugin, Molecular mechanisms of lung cell activation induced by cyclic stretch, Crit. Care Med. 31 (4) (2003) S200–S206. [8] E. Denny, R.C. Schroter, A model of non-uniform lung parenchyma distortion, J. Biomech. 39 (4) (2006) 652–663. [9] M. Felici, et al., Diffusional screening in real 3D human acini—a theoretical study, Respir. Physiol. Neurobiol. 145 (2–3) (2005) 279–293. [10] L. Harrington, G.K. Prisk, C. Darquenne, Importance of the bifurcation zone and branch orientation in simulated aerosol deposition in the alveolar zone of the human lung, J. Aerosol Sci. 37 (1) (2006) 37–62. [11] H. Kitaoka, et al., A 4-dimensional model of the alveolar structure, J. Physiol. Sci. 57 (3) (2007) 175–185. [12] B.S. Ma, C. Darquenne, Aerosol deposition characteristics in distal acinar airways under cyclic breathing conditions, J. Appl. Physiol. 110 (5) (2011) 1271–1282. [13] B. Suki, et al., In silico modeling of interstitial lung mechanics: implications for disease development and repair, Drug Discov. Today:Dis. Models 4 (3) (2007) 139–145. [14] A.J. Swan, M.H. Tawhai, Evidence for minimal oxygen heterogeneity in the healthy human pulmonary acinus, J. Appl. Physiol. 110 (2) (2011) 528–537. [15] J. Sznitman, et al., Respiratory flow phenomena and gravitational deposition in a three-dimensional space-filling model of the pulmonary acinar tree, J. Biomech. Eng.-Trans. ASME 131 (2009) 3. [16] S.M. Rausch, et al., Local strain distribution in real three-dimensional alveolar geometries, Annu. Biomed. Eng. 39 (11) (2011) 2835–2843. [17] A. Tsuda, et al., Finite element 3D reconstruction of the pulmonary acinus imaged by synchrotron X-ray tomography, J. Appl. Physiol. 105 (3) (2008) 964–976. [18] E. Denny, R.C. Schroter, A mathematical model for the morphology of the pulmonary acinus, J. Biomech. Eng. 118 (2) (1996) 210–215. [19] Y.C. Fung, A model of the lung structure and its validation, J. Appl. Physiol. 64 (5) (1988) 2132–2141.
[20] H. Kumar, et al., The effects of geometry on airflow in the acinar region of the human lung, J. Biomech. 42 (11) (2009) 1635–1642. [21] J. Sznitman, Respiratory microflows in the pulmonary acinus, J. Biomech. 46 (2) (2013) 284–298. [22] H. Parameswaran, A. Majumdar, B. Suki, Linking microscopic spatial patterns of tissue destruction in emphysema to macroscopic decline in stiffness using a 3D computational model, PLOS Comput. Biol. 7 (2011) 4. [23] J.B. West, Regional differences in the lung, Chest 74 (4) (1978) 426–437. [24] E. Roan, C.M. Waters, What do we know about mechanical strain in lung alveoli? Am. J. Physiol.-Lung Cell. Mol. Physiol. 301 (5) (2011) L625–L635. [25] B. Suki, D. Stamenovic, R. Hubmayr, Lung parenchymal mechanics, Compr. Physiol. 1 (3) (2011) 1317–1351. [26] D. Haberthur, et al., Visualization and stereological characterization of individual rat lung acini by high-resolution X-ray tomographic microscopy, J. Appl. Physiol. 115 (9) (2013) 1379–1387. [27] T. Sera, et al., Murine pulmonary acinar mechanics during quasi-static inflation using syncrotron refraction-enhanced CT, J. Appl. Physiol. 115 (2013) 219–228. [28] D.M. Vasilescu, et al., Assessment of morphometry of pulmonary acini in mouse lungs by nondestructive imaging using multiscale microcomputed tomography, Proc. Natl. Acad. Sci. USA 109 (42) (2012) 17105–17110. [29] K.S. Burrowes, et al., Towards a virtual lung: multi-scale, multi-physics modelling of the pulmonary system, Philos. Trans. R. Soc. A-Math. Phys. Eng. Sci. 366 (1879) (2008) 3247–3263. [30] H. Kitaoka, S. Tamura, R. Takaki, A three-dimensional model of the human pulmonary acinus, J. Appl. Physiol. 88 (6) (2000) 2260–2268. [31] A. Okabe, et al., Spatial Tessellations: Concepts and Applications of Voronoi Diagrams. Wiley Series in Probability and Statistics, Wiley, England, 2000. [32] D.S. Johnson, et al., Optimization by simulated annealing—an experimental evaluation.1. graph partitioning, Oper. Res. 37 (6) (1989) 865–892. [33] S. Kirkpatrick, C.D. Gelatt, M.P. Vecchi, Optimization by simulated annealing, Science 220 (4598) (1983) 671–680. [34] A. Leeming, R. Schroter, A model morphology of the pulmonary acinus, Proc. Inst. Mech. Eng. Part H-J. Eng. Med. 222 (H4) (2008) 429–437. [35] C.R. Ethier, C.A. Simmons, Introductory Biomechanics From Cell to Oranisms, Cambridge University Press, New Yourk, 2007. [36] C.B. Barber, D.P. Dobkin, H. Huhdanpaa, The Quickhull algorithm for convex hulls, ACM Trans. Math. Softw. 22 (4) (1996) 469–483. [37] L. Xiao, et al., A semiautomatic segmentation algorithm for extracting the complete structure of acini from synchrotron micro-CT images, Comput. Math. Methods Med. 575086 (2013) 1–10. [38] J.E. Hansen, E.P. Ampaya, Human air space shapes, sizes, areas, and volumes, J. Appl. Physiol. 38 (6) (1975) 990–995. [39] J.E. Hansen, et al., Branching pattern of airways and air spaces of a single human terminal bronchiole, J. Appl. Physiol. 38 (6) (1975) 983–989. [40] R.R. Mercer, J.M. Laco, J.D. Crapo, 3-Dimensional reconstruction of alveoli in the rat lung for pressure-volume relationships, J. Appl. Physiol. 62 (4) (1987) 1480–1487. [41] R.R. Mercer, M.L. Russell, J.D. Crapo, Alveolar septal structure in different species, J. Appl. Physiol. 77 (3) (1994) 1060–1066. [42] H. Parameswaran, et al., Three-dimensional measurement of alveolar airspace volumes in normal and emphysematous lungs using micro-CT, J. Appl. Physiol. 107 (2) (2009) 583–592.
©2015 Elsevier
Contents lists available at ScienceDirect
Computers in Biology and Medicine journal homepage: www.elsevier.com/locate/cbm
Mathematical model of a heterogeneous pulmonary acinus structure Kenichiro Koshiyama n, Shigeo Wada Graduate School of Engineering Science, Osaka University, Japan
art ic l e i nf o
a b s t r a c t
Article history: Received 30 October 2014 Accepted 31 March 2015
The pulmonary acinus is a gas exchange unit distal to the terminal bronchioles. A model of its structure is important for the computational investigation of mechanical phenomena at the acinus level. We propose a mathematical model of a heterogeneous acinus structure composed of alveoli of irregular sizes, shapes, and locations. The alveoli coalesce into an intricately branched ductal tree, which meets the space-filling requirement of the acinus structure. Our model uses Voronoi tessellation to generate an assemblage of the alveolar or ductal airspace, and Delaunay tessellation and simulated annealing for the ductal tree structure. The modeling condition is based on average acinar and alveolar volume characteristics from published experimental information. By applying this modeling technique to the acinus of healthy mature rats, we demonstrate that the proposed acinus structure model reproduces the available experimental information. In the model, the shape and size of alveoli and the length, generation, tortuosity, and branching angle of the ductal paths are distributed in several ranges. This approach provides a platform for investigating the heterogeneous nature of the acinus structure and its relationship with mechanical phenomena at the acinus level. & 2015 Elsevier Ltd. All rights reserved.
Keywords: Pulmonary biomechanics Biological microstructure Pulmonary physiology Gas exchange Computational geometry Structure modeling Rat acinus
1. Introduction The pulmonary acinus is a gas exchange unit that is located distal to the terminal bronchioles and consists of a complex microstructure of ducts and alveoli [1–3]. The mechanical phenomena at the acinus level, i.e., acinar flow, gas diffusion, particle depositions, or mechanical responses of lung parenchyma to stretches, underlie the physiological behavior of the lungs. A detailed understanding of the pulmonary acinus is required for recent developments in medical engineering techniques, e.g., pulmonary drug delivery [4,5] or mechanical ventilation protocols [6,7]. To date, pulmonary researchers have investigated acinus phenomena using theoretical and numerical approaches [8–15]. The structure model of the acinus (i.e., the shape and size of the alveoli and the patterns of the ductal tree) affects and restricts the numerical results, and hence the choice of structure model is important. Various image-based [16,17] and mathematical models [11,18–21] have been proposed. Image-based models represent part of the real acinus structure, but are sensitive to the image resolution and segmentation algorithm. Hence, their application to numerical analysis is still restricted. In contrast, mathematical models can easily mimic general characteristics of the acinus
n Correspondence to: Graduate School of Engineering Science, Osaka University, Toyonaka, 560-8531, Japan. Tel.: þ 81 6 6850 6174; fax: þ81 6 6850 6173. E-mail address: [email protected] (K. Koshiyama).
http://dx.doi.org/10.1016/j.compbiomed.2015.03.032 0010-4825/& 2015 Elsevier Ltd. All rights reserved.
structure. Two types of models are frequently used for further numerical analysis. The first type is composed of straight tubes that enable regular sphere caps or polyhedra to coalesce on the side wall (i.e., duct-based model) and is typically used to analyze acinar flows [21]. The second type of model is formed from polyhedra of regular size and shape (i.e., polyhedra-based model) and is used to analyze the mechanical response of lung parenchyma to stretching [22]. Conventional models rely on the average morphological characteristics of the acinus structure. Consequently, they are too basic to capture the heterogeneous nature of the alveolar and ductal patterns observed in real acinus structures. Computations using simplified acinus models provide approximated features of mechanical phenomena, but allow only a limited investigation of phenomena related to the heterogeneous nature of the acinus structure. In fact, variations of mechanical phenomena in the lungs (i.e., regional differences [23]) have long been recognized. Advances in mechanobiology mean that the detailed resolution of the strain field in alveolar walls is of great interest to researchers in lung cell mechanics [24]. Considering lung parenchyma mechanics, the pathway for stress transmission from organ to cells or cytoskeletons and vice versa is through the extracellular matrix. The parenchyma architecture i.e., detailed acinus structure, is one of important factors affecting such mechanotransduction [13,25]. The heterogeneity of the acinus structure on acinar flow remains important in biofluidics [21]. Furthermore, the acinus structure segmented from high-resolution CT images has been analyzed to micrometer resolution [26–28], which has the potential to provide
26
K. Koshiyama, S. Wada / Computers in Biology and Medicine 62 (2015) 25–32
quantitative data about the heterogeneity of the acinus structure. Thus, to understand the heterogeneous nature of the acinus structure and its relationship to mechanical phenomena in the lungs, it is vital to develop a parametric heterogeneous acinus structure model and investigate the potential and limitations of the modeling technique. In this paper, we propose a method for generating a heterogeneous acinus structure model composed of alveoli of irregular sizes, shapes, and locations that coalesce into an intricately branched ductal tree. Instead of including all the details of the actual acinus structure, we focus on the gas exchange unit of the 3D space-filling structure, represented by an assemblage of polyhedra. To our knowledge, several groups have developed heterogeneous acinus structure models [29,30], but they have focused on the heterogeneity of the ductal tree. A distinctive feature of our model is that we parametrically and explicitly include the distributions of both alveolar and ductal structures, without significantly departing from the typical polyhedra-based acinus model [18,19]. This is achieved by introducing Voronoi and Delaunay tessellations [31] and simulated annealing optimization [32,33]. The abundance of quantitative and qualitative morphological data about the acinus structure [2,34] allow us to verify our approach by modeling the acinus structure of a rat and comparing the model output with previous results. Our paper is organized as follows. In Section 2, we describe the basic methodology of heterogeneous modeling, present the model setting for the rat acinus structure, and introduce the structural parameters analyzed. Section 3 describes the sensitivity of the average structural parameters to the modeling method, distributions and variations of alveolar and ductal structures, segmentation of the 3D structure, model predictions, and limitations. Finally, Section 4 summarizes the paper.
2. Methods The approach presented here is a unified version of the polyhedra-based acinus structure modeling first presented by Denny and Schoroter [18]. The polyhedra-based modeling regards each polyhedron as an airspace region and the faces as tissue regions, i.e., alveolar or ductal walls. The approach is composed of two modeling steps: the first generates an assemblage of polyhedra, and the second generates a ductal tree by eliminating the common walls between these polyhedra based on the optimized ductal tree when the gas exchange efficiency in the acinus is considered as the cost for the optimization. In the following, we describe the extension of this approach to express the heterogeneity of the acinus structure. We illustrate the 2D case for simplicity (Fig. 1), although the method can be extended to three dimensions.
that when the seed coordinates are set to the body-centered cubic lattice in a 3D space, the structure of the resultant Voronoi polyhedron becomes a truncated octahedron, which is a typical shape for alveolar or ductal airspaces [18,19,29]. 2.2. An intricately branched ductal tree A tree structure is generated by connecting the neighboring seeds. To consistently obtain the neighboring seeds in the heterogeneous system, we introduce the Delaunay tessellation (Fig. 1C) from the mathematical background of the Voronoi tessellation [31]. To generate a fully isolated structure, we do not consider connections to seeds in polygons whose edges open outward (boundary polygons) (Fig. 1B). A seed near the center of the top edge of a given region is designated as the entrance seed, and subsequent connections are randomly selected from the neighboring seeds (Fig. 1D). The tree generation obeys two geometrical constraints based on experimental findings from an acinar ductal tree [2,18]: (i) each seed has a unique path to the entrance seed (i. e., no cycle); and (ii) the number of daughter branches from a seed is less than 3, whereas the number of branches to the terminal seed, i.e., to the alveoli, is unlimited. Because the acinus is a gas exchange unit, we control the geometrical components related to gas exchange efficiency by optimizing the tree structure based on the following criteria. The seeds are numbered, and a solution in the optimization provides each seed with information about one parent (upstream) seed number and the daughter (downstream) seed numbers (i.e., connective information). The neighbor of a solution is obtained by randomly changing the parent seed of a randomlyselected seed. As the gas exchange in real acini involves a complex set of factors (e.g., gas diffusion and convection in airways, gas exchange between airways and surrounding lung capillary networks, blood circulation or respiration conditions [35]), the selection of cost for the optimization is arbitrary. We assume that efficient gas exchange in acini is related to a shortening of ductal path lengths for gas diffusion and an increase in the number of alveoli for the region of the gas exchange which can be calculated from connective information with less computational costs. The same criterion has been used previously [18,34]. Accordingly, the cost to be minimized during optimization includes the cost C 1 of maximizing the number of alveoli (i.e., to minimize the number of ducts) C 1 ¼ Nt Na ;
where N t is the total number of seeds and N a is the total number of terminal seeds, as calculated from the connective information; and the cost C 2 of minimizing the average path length from the entrance to each terminal seed
2.1. An assemblage of airspace regions of irregular shapes and sizes C2 ¼ Firstly, seeds of polygons (polyhedra in 3D) are scattered in a given (square) space (Fig. 1(A)). The spatial coordinates are given by sets of pseudo-random real numbers taken from a uniform distribution between 0 to the distance in one direction of the space. We set the number density of seeds and the minimum distance as the conditions for seed generation. Computationally, for a given seed density, decreasing the minimum distance increases free space to generate seed coordinates, hence allowing various seed distributions. Polygons filling a given space are generated based on seeds by the digitized Voronoi tessellation technique, which divides a region by allocating points to the closest seeds [31] (Fig. 1A and B). Accordingly, each seed has information about edges of its own polygon (faces of its own polyhedron in 3D). It should be noted
ð1Þ
Na 1 X l; Na i ¼ 1 i
ð2Þ
where li is the path length to the ith terminal seed, which is the sum of the lengths between the connected seeds' coordinates in the system. The total cost C is defined as C ¼ ω C 1' þ ð1 ωÞC 2' ;
ð3Þ
where ω is the weight of each cost component, 0 r ω r 1, and the apostrophe denotes that each cost is normalized according to its values at ω ¼ 0 or 1. For the optimization, we introduce a simulated annealing algorithm, where the optimized tree structure is searched by modifying the connections between seeds and comparing the costs before and after the modification. Specifically, the parent of
K. Koshiyama, S. Wada / Computers in Biology and Medicine 62 (2015) 25–32
27
Fig. 1. Schematic diagram of heterogeneous acinus structure modeling: random seed distribution (A), Voronoi diagram (B), Delaunay diagram (C), initial ductal tree structure (D), optimized ductal tree structure (E), and resulting acinus structure model (airspace is shown in gray) (F).
Table 1 Rat acinus modeling parameters and their reproducibility. Dmin [μm]
Lx [mm]
ρ Vt F [-/mm3] [106 μm3] [%]
Va [μm3]
RV t [%] RF [%] RVa [%]
50 40 30 Experiment [2]
0.97 0.98 0.99 –
5478 5312 5153 –
196 190 187 194
99 99 98 100
527 527 532 522
76 76 75 80
95 94 94 100
99 98 98 100
Dmin : minimum distance between seeds; Lx : distance in one direction of the cubic system; ρ: number density of seeds; V t : average total airspace volume; F: average fraction of airspace; and V a : average alveolar volume. RVt , RF , and RV a denote the reproducibility of V t , F, and V a (see Section 2.4).
a randomly chosen seed is changed to a neighboring seed on the Delaunay sides, and then the new cost is calculated via Eq. (3) with the new tree structure. When the new cost is smaller than the previous one, the new tree structure is accepted; even when the new cost is larger than the previous one, it may be accepted on statistical terms (i.e., under the Metropolis criterion). The probability of acceptance of a larger cost decreases with the number of modifications, and the modification is iterated until all costs have converged. After the optimized tree structure was obtained (Fig. 1E), the polygonal edges pierced by the optimized tree are eliminated. As each seed has information about edges of its own polygon, eliminating the common edges of connecting seeds results in all the polygons in the system being vented from the entrance (Fig. 1F). 2.3. 3D acinus model settings To visualize and analyze the 3D geometry of the acinus model, a volume rendering technique is required. In practice, this is achieved by forming the 3D digitized Voronoi polyhedra using
an object of arbitrary shape composed of a number of fine cubic building blocks (voxels). The dimensions of the object, voxel size, number of voxels representing the acinus tissue (i.e., alveolar wall thickness), and number of seeds are related to the acinar and alveolar sizes. These parameters are set so that the total airspace volume, average alveolar volume, and volume density of the airspace in the model match those determined experimentally. The three volume-related data are the modeling condition for 3D acinus modeling. For simplicity, we assume that the acinus is cubic, all alveoli belong to a single acinus, and the number of seeds, the voxel size, and the wall thickness are constant. Furthermore, the seed coordinates are given by a pseudo-random number generator with some minimum distance between seeds, by which we attempt to control spatial distributions of seeds. For the rat acinus structure, the number of seeds, the voxel size, and wall thickness are set to 5000, 2 μm, and 4 μm, respectively; other modeling parameters are listed in Table 1. For generating irregular seed distribution, sets of pseudo-random numbers, taken from a uniform distribution in the range ð0; 1Þ, are scaled by the distance in one direction of the cubic system and are used as spatial coordinates of seeds. The minimum distances are 30, 40, and 50 μm. Note that, in practice, the modeling parameters for each minimum distance are set empirically to fulfill the volumerelated data from experiments. Accordingly, the resultant number densities of seeds for different minimum distances differ slightly (Table 1), but controlling seed distributions with minimum distance work with the given conditions and assumptions (see Section 3.2). #Although other structural information is available for the rat acinus (e.g., number of alveoli per acinus, average path length, etc. [2]), volume-related information is relatively reliable, as it can be obtained by simply counting the voxels in acinus images, whereas detailed alveolar or ductal structural information requires interpretation (e.g., skeletonization for the average path length, and estimation from small airway samples or 2D images for the
28
K. Koshiyama, S. Wada / Computers in Biology and Medicine 62 (2015) 25–32
structures increase with the weight, as expected, and exhibit a plateau region from 0:2 to 0:8, which corresponds to the previous regular polyhedra-based model (Fig. 4 in [34]). Additionally, the trend is insensitive to the initial seed distribution and the minimum distance between seeds (Fig. S2 in Supplementary information file A). The maximum number of generations, which has been shown to be important in verifying the path structure [34], and other path structures also exhibit a plateau region when the weight is not close to 0 or 1 (Fig. 2B and Fig. S3 in Supplementary information file A). These results suggest that, as long as we consider the two costs in Eq. (3) for gas exchange efficiency, the average path structure is essentially invariant. This conclusion is confirmed by the statistics of experimental results [2] and the previous regular structure model [34]. Although this insensitivity to weight (away from the extreme values) is evident in the current and previous [34] polyhedrabased models, there are slight differences in the invariant values. The current model gives approximately 2800 alveoli in a maximum of 14 generations, whereas the previous model suggests around 2200 alveoli in 12 generations. This difference can be explained by the explicit modeling of the wall thickness and alveolar size distributions in the current model, as the total number of polyhedra ( 3200 after eliminating boundary seeds) for the volume-related modeling condition (Table 1) is larger than that of the previous model ( 2700) [34].
number [27]). Therefore, in the current study, we use volumerelated data as the modeling conditions. For the simulated annealing, we introduce the Johnson–Aragon–McGeoch–Schevon algorithm [32] (detailed algorithm and parameters are given in Supplementary information file A). To verify the effects of the weight in Eq. (3) on the resulting acinus structure, we change the weight from 0 to 1. For each weight, we analyze the data from three different initial seed distributions. The modeling and analysis use in-house code written in C and Fortran 90 with the aid of the Qhull computational geometry algorithms [36] and the Mersenne Twister random seed generator. The source code is available from the authors upon request. To visualize the 3D geometry, we render sets of sliced images and reconstruct them using the general visualization tool Amillas. 2.4. Analysis Based on the definition given in previous studies [1,34], we regard the path from a branching seed to the next branching seed as one generation, ignoring the path to the terminal seed (alveoli). The complexity of the ductal tree is characterized by the distributions of the tortuosity and of the branching angle. The tortuosity is estimated as the ratio between the actual path length and the linear distance from the entrance seed to the terminal branching seed, where the actual path length is the sum of the path lengths per generation. The branching angle is that between the direction to the periphery of the parent's path and that of the daughter's path. The reproducibility R of the experimental information is calculated as R ¼ 1 jV V vj, where v is the average value given by the model and V is the experimental value. R is a distance measure, i.e., absolute deviation, between the model value and the experimental value, normalized by the experimental value. In order to mean that, at R ¼ 1:0(100%), the model value is equivalent to the experimental value, the distance is shifted one unit to the positive. R is only valid for 0 r v r 2Vas R becomes negative for v 4 2V.
3.2. Structural distributions for ω ¼ 0:5 and effects of minimum distance Because of the insensitivity of the average structure to the cost weighting, it is worth investigating the distributions of various structural parameters for a particular value (ω ¼ 0:5). One important characteristic of the polyhedra-based acinus model is that it expresses certain path structure distributions, i.e., path length, number of generations, tortuosity, and branching angles [18] (see Fig. 3C–F). In the current model, we extend this feature to express the distributions of alveolar shapes and sizes (Fig. 3A and B). By changing the initial seed distribution, we can create structure models in which the shape of each distribution is different while the average values remain almost constant. Fig. 4 shows the effects of the minimum distance Dmin on the alveolar shape and volume. For all Dmin , the peak value of the number of polyhedra faces is 14, which corresponds to the number of faces in the truncated octahedron often used to give a mathematical representation of the airspace in acinus [19]. The distributions of the number of faces and volume of alveoli become narrower as the minimum distance increases. Interestingly, the path structure distributions are insensitive to changes in this
3. Results and discussion 3.1. Effect of weight on the acinus structure To verify the sensitivity of the optimized structure to the weighting parameter in the heterogeneous acinus modeling, we first examined the sensitivity of the number of alveoli, average path length, and maximum number of generations. Fig. 2A shows the effect of ω on the number of alveoli and the average path length, which are directly optimized via costs C 1 and C 2 . Both
2500
3 No. of alveoli Av. path length
2000
2
1500
1
1000
0
0.2
0.4
0.6
ω[−]
0.8
1
0
Number of generations [ - ]
4
Average path length [ mm ]
Number of alveoli [ - ]
20 3000
18 16 14 12 10 0
0.2
0.4
0.6
0.8
1
ω[−]
Fig. 2. Effect of cost weight on the number of alveoli and average path length (A) and the maximum number of generations (B) for Dmin ¼ 40. Average values from three different seed distributions are shown. Error bars in Fig. 2B denote the standard deviation.
K. Koshiyama, S. Wada / Computers in Biology and Medicine 62 (2015) 25–32
30
30 20
30 20
Model 2
10
Model 1
10
Frequency [%]
Frequency [%]
20
Model 1
10
30 20
30 20
Model 2
10 30 20
Model 3
10 0
0
5
10
15
20
25
Model 3
10 0
30
0
3.0
4.0 2 ×10
30
20
20
Model 1
Model 1
10
Frequency [%]
10
Frequency [%]
2.0
Alveolar volume [ µm ]
30
30 20
Model 2
10 30 20
30 20
Model 2
10 30 20
Model 3
10
Model 3
10 0
0.5
1.0
1.5
0
2.0
0
Path length [ mm ]
10
15
20
30
20
20
Model 1
30 20
Model 2
10
Model 1
10
Frequency [%]
10
30
30 20
Model 2
10 30 20
20
Model 3
10 0 1.0
5
Number of generations [ - ]
30
Frequency [%]
1.0
3
Number of faces [ - ]
0
29
1.5
2.0
2.5
Model 3
10 3.0
Tortuosity [ - ]
0
0
45
90
135
180
Angle [ - ]
Fig. 3. Representative results about distributions of the number of faces (A), alveolar volume (B), path length (C), number of generations (D), tortuosity (E), and path angle (F) for Dmin ¼ 40 with ω ¼ 0:5. Model results generated from three different initial seed distributions are labeled as Model 1, 2, and 3.
minimum distance (Fig. S4 in Supplementary information file A). Because we have used a cubic acinus model with almost the same total volume for different Dmin (Table 1), we suspect that the shape of the acinus might affect the path structure distributions. We are now developing an ellipsoidal acinus model, and are analyzing the shape and volume information of individual acini from CT images [37]. The effects of acinar shape on the path structure will be investigated in more detail in the near future. 3.3. Variations of acinus structure during optimization The acinus models presented in Sections 3.1 and 3.2 and in previous studies [18,34] are based on the optimized tree structure. During optimization, the tree experiences various structures.
Fig. 5A presents representative results for the variation in average path length and number of alveoli. The two structural parameters change on each iteration to fulfill the modeling conditions (Table 1). Depending on the weight and initial seed distribution, the paths of the variations change (Fig. S5). For the modeling conditions of volume-related experimental data (Table 1), the average path length and the number of nonterminal seeds are minimized regardless of experimental values. As experimental values for the number of alveoli and average path length are available ( 1542 and 1.287 mm, respectively [2]), the reproducibility of these parameters can be examined along the path (Fig. 5B). Interestingly, the reproducibility approaches 1, before moving away during optimization. This indicates that, in the given condition, the actual rat acinus structure is not fully
30
K. Koshiyama, S. Wada / Computers in Biology and Medicine 62 (2015) 25–32
30
40 30
20 Dmin = 30 mm
Frequency [ % ]
Frequency [ % ]
10 30 20 Dmin = 40 mm 10 30 20
0
0
5
15
10
20
25
30 Dmin = 40 mm
20 10 40 30 20 10
Dmin = 50 mm
10
Dmin = 30 mm
20 10 40
30
0
Dmin = 50 mm
0
1.0
2.0
3.0
4.0
Alveolar volume [ µm 3 ]
Number of faces [ - ]
5.0
×10 2
Reproducibility
Average path length [ mm ]
Fig. 4. Effects of the minimum distance between seeds on the distributions of the number of faces (A) and alveolar volume (B). Distributions were calculated using data obtained from three models for each Dmin under ω ¼ 0:5.
6 5 4 3 2 1 0 1000
the distributions of alveolar and ductal structural parameters (Figs. 2–5) mean that the current model provides a platform for computational investigations of the heterogeneous nature of the acinus structure and its relationships with mechanical phenomena at the acinus level.
Initial tree Optimized tree
1500
2000 Number of alveoli [ - ]
2500
3000
1.0 Average path length
0.8 0.6 0.4 0.2 0
Number of alveoli Total
10 5 10 6 Number of iterations [ - ]
107
Fig. 5. Representative relationship between average path length and number of alveoli during optimization for ω ¼ 0:6 (A), and the change in reproducibility (B). In (B), the black solid line represents the reproducibility of the number of alveoli, gray dashed line represents average path length, and black bold line denotes total reproducibility.
optimized for gas exchange efficiency, at least with the given cost functions and the modeling assumptions in Section 2, and other mechanical factors may be involved in the morphology.
3.4. Segmentation of heterogeneous acinus structure models Segmentations of the airspace or tissue region are valuable for further numerical calculations. For this purpose, we explicitly modeled the alveolar wall and performed an image rendering (Fig. 6). After determining the ductal tree, and eliminating the walls that the ductal tree pierces, in a 3D cube (Fig. 6A and B), slices of the binarized images (in which tissue voxels and airspace voxels are shown in black and white, respectively) are stacked (Fig. 6C). This is a similar approach to the generation of 3D models from 2D images [2,28]. When the black voxels are reconstructed by general image processing tools (Fig. 6D), the model can be readily rendered as a 3D mesh for mechanical analyses, e.g., the mechanical response of the parenchyma to stretches, and when the white voxels are reconstructed and meshed (Fig. 6E), the model can be used to analyze gas diffusion or acinar flow. Part of the acinus (alveolar ducts or alveolar sac) can also be reconstructed (Fig. 6F). The variety of the average structure (Fig. S5) and
3.5. Model predictions and limitations A distinctive feature of our model is that we parametrically and explicitly include the distributions of both alveolar and ductal structures (Figs. 3–6). In addition to the comparison with the average data (Table 1 and Fig. 5B), we compare the distributions (Figs. 3 and 4) with available experimental data. Hansen and Ampaya [38] analyzed airspace regions of human acini and presented six shapes for alveolar airspace, assuming symmetrical shapes with curved surfaces. Fung [19] proposed 14hedron is a simple alveolar shape. Although we did not assume symmetrical shapes, the peak value of the number of polyhedral faces was 14 and they were distributed from 8 to 25, depending on the minimum distance (Fig. 4A). Mercer and Crapo [2] reported that the alveolar volume is 194 721 μm3 (average 7 standard error) for similar acinar conditions used for our modeling (Table 1). Although they provided only the standard error, the standard deviation might be large, considering the number of alveoli they analyzed (sixty alveoli from different lungs). The distributions of alveolar volume and shape can be controlled in our modeling by changing the minimum distance between seeds (Fig. 4), while a regular model provides only airspace region of regular volume and shape [18]. Among the pathway structure distributions, we verify the length distribution (Fig. 3C) because the frequency distributions for four individual acini have been reported by Mercer and Crapo [2]. From their experiment, regardless of acini analyzed, the lengths range from 0.2 to 2.5 mm and the maximum relative frequency is under 30%. The similar range of the path length for average acinus has been presented by Rodriguez, et al. [3]. Likewise, the lengths from our model are within the range of the above experimental data (Fig. 3C). In addition, the experimental data about different acini suggest that the shape of the distribution depends on individual acini [2]. The difference of the distribution shape can be obtained in our models for different seed distributions (Fig. 3C). The variety of alveolar and pathway structures in acinus has been widely observed in experiments [1–3,27,38–42]. The average data and their distributions depend on, species [41], inflation level [27,40], region in lungs [23], pathological conditions [42], or ages
K. Koshiyama, S. Wada / Computers in Biology and Medicine 62 (2015) 25–32
31
Fig. 6. Representative 3D heterogeneous acinus rendering: seeds distributed in a cube (A), optimized path tree (B), stack of 2D images (C), wall region rendering (D), airspace rendering (E), and part of the acinus model (F). See the text for more details.
[28]. Strictly, methodological issues (fixation of samples, imaging devices, segmentation algorithm, etc.) may also affect the data because of the delicate microstructure of acinus. Hence, great care is needed for the quantitative comparison of model predictions with experimental data. Nevertheless, as shown in Figs. 3–6, the architecture of acinus model presented here includes distributions of, not only pathway structures, but also alveolar shape, size, and locations, while the regular model has the distribution of pathway structures and alveolar location only [18]. Furthermore, we can parametrically control the distributions alveolar shape and size by setting minimum distance between seeds (Fig. 4). Besides the advantage, our approach has some limitations. First, the duct size may vary with the generation of ducts in the actual acinus [2,3], whereas our model assumes the size is independent of generation. Second, the alveolar shape may be restricted in the actual acinus [38], partly due to mechanical balance in acinus, whereas it is unlimited in our model. The first limitation is common in polyhedra-based models, and restricts our modeling to acini that do not include apparent respiratory bronchioles. Applied to human acinus, our model may correspond to the 1/8 subacinus [1], and hence, volume-related data should be taken from this region. The second point may affect the reproducibility (Fig. 5B). Because our model is a static model and does not involve any mechanics, further improvements should be made to e.g., the cost functions used in the current model. Alternatively, it may be interesting to investigate the effects of mechanics of alveolar walls (stiffness to the septal walls, negative pressure boundary, surface tension effects, etc. [25]) on the shape and size distributions by mechanical analysis using our structure models [8,16,25]. The assumptions used in the modeling, e.g., cubic acinar shape, constant wall thickness and number of seeds for setting parameters, etc., are also the limitations and should be addressed in the future studies. The resolution of in situ CT acinus imaging has reached the scale of the alveolar wall in various inflation levels [17,27], and our group has developed a semiautomatic segmentation algorithm that has the potential to provide the morphology of individual acini and a cluster of acini stemming from the same conducting airway [37]. Even in such detailed images, identifying ducts or alveoli is still challenging, but more sophisticated statistical data about alveolar and ductal distributions of individual acini will facilitate and complement our modeling, and vice versa.
4. Summary We have demonstrated that the combined use of Voronoi and Delaunay tessellations and the optimization of the ductal tree with simulated annealing can parametrically model a heterogeneous acinus structure composed of alveoli of irregular sizes, shapes, and locations that coalesce into an intricately branched ductal tree. The modeling condition considered average acinar and alveolar volume characteristics from published experimental information on rat acinus. In the model, the shape and size of alveoli and the length, generation, tortuosity, and branching angle of the ductal paths were distributed in several ranges. The distribution of the shape and size of the alveoli depends on the minimum distance between airspace regions in acinus models, while those of the path structure were found to be insensitive in cubic acinus models. The optimized tree structure did not fully reproduce the experimental values of the number of alveoli and average path length, indicating that the actual rat acinus structure is not fully optimized for gas exchange efficiency and that other mechanical factor may be involved in the morphology. The proposed technique can model the tissue and airspace of acinus without significantly departing from the typical acinus model used for various mechanical analyses. Therefore, the heterogeneous acinus model developed here can be used to investigate the heterogeneous nature of the acinus structure and its relationships with mechanical phenomena at the acinus level using computational approaches.
Conflict of interest statement The authors have no conflict of interests to declare with regard to this study.
Acknowledgments We would especially like to thank Prof. R.C. Schroter and Dr. T. Sera for their useful comments and suggestions. We also thank M. Yasui, A. Mizuno, T. Nishimura, and L. Xiao for their technical assistance. This work was partly supported by the Japanese Ministry of Education, Culture, Sports, Science and Technology (24240073).
32
K. Koshiyama, S. Wada / Computers in Biology and Medicine 62 (2015) 25–32
Appendix A. Supporting information Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j.compbiomed.2015. 03.032. References [1] B. Haefeli-Bleuer, E.R. Weibel, Morphometry of the human pulmonary acinus, Anat. Rec. 220 (4) (1988) 401–414. [2] R.R. Mercer, J.D. Crapo, Three-dimensional reconstruction of the rat acinus, J. Appl. Physiol. 63 (2) (1987) 785–794. [3] M. Rodriguez, et al., Pulmonary acinus—geometry and morphometry of the peripheral airway system in rat and rabbit, Am. J. Anat. 180 (2) (1987) 143–155. [4] N.R. Labiris, M.B. Dolovich, Pulmonary drug delivery. Part I: physiological factors affecting therapeutic effectiveness of aerosolized medications, Br. J. Clin. Pharmacol. 56 (6) (2003) 588–599. [5] J.S. Patton, J. Bukar, S. Nagarajan, Inhaled insulin, Adv. Drug Deliv. Rev. 35 (2–3) (1999) 235–247. [6] K.W. Delank, et al., Diagnosis and treatment modalities of sinunasal inverted papillomas, Laryngo-Rhino-Otol. 79 (4) (2000) 226–232. [7] J. Pugin, Molecular mechanisms of lung cell activation induced by cyclic stretch, Crit. Care Med. 31 (4) (2003) S200–S206. [8] E. Denny, R.C. Schroter, A model of non-uniform lung parenchyma distortion, J. Biomech. 39 (4) (2006) 652–663. [9] M. Felici, et al., Diffusional screening in real 3D human acini—a theoretical study, Respir. Physiol. Neurobiol. 145 (2–3) (2005) 279–293. [10] L. Harrington, G.K. Prisk, C. Darquenne, Importance of the bifurcation zone and branch orientation in simulated aerosol deposition in the alveolar zone of the human lung, J. Aerosol Sci. 37 (1) (2006) 37–62. [11] H. Kitaoka, et al., A 4-dimensional model of the alveolar structure, J. Physiol. Sci. 57 (3) (2007) 175–185. [12] B.S. Ma, C. Darquenne, Aerosol deposition characteristics in distal acinar airways under cyclic breathing conditions, J. Appl. Physiol. 110 (5) (2011) 1271–1282. [13] B. Suki, et al., In silico modeling of interstitial lung mechanics: implications for disease development and repair, Drug Discov. Today:Dis. Models 4 (3) (2007) 139–145. [14] A.J. Swan, M.H. Tawhai, Evidence for minimal oxygen heterogeneity in the healthy human pulmonary acinus, J. Appl. Physiol. 110 (2) (2011) 528–537. [15] J. Sznitman, et al., Respiratory flow phenomena and gravitational deposition in a three-dimensional space-filling model of the pulmonary acinar tree, J. Biomech. Eng.-Trans. ASME 131 (2009) 3. [16] S.M. Rausch, et al., Local strain distribution in real three-dimensional alveolar geometries, Annu. Biomed. Eng. 39 (11) (2011) 2835–2843. [17] A. Tsuda, et al., Finite element 3D reconstruction of the pulmonary acinus imaged by synchrotron X-ray tomography, J. Appl. Physiol. 105 (3) (2008) 964–976. [18] E. Denny, R.C. Schroter, A mathematical model for the morphology of the pulmonary acinus, J. Biomech. Eng. 118 (2) (1996) 210–215. [19] Y.C. Fung, A model of the lung structure and its validation, J. Appl. Physiol. 64 (5) (1988) 2132–2141.
[20] H. Kumar, et al., The effects of geometry on airflow in the acinar region of the human lung, J. Biomech. 42 (11) (2009) 1635–1642. [21] J. Sznitman, Respiratory microflows in the pulmonary acinus, J. Biomech. 46 (2) (2013) 284–298. [22] H. Parameswaran, A. Majumdar, B. Suki, Linking microscopic spatial patterns of tissue destruction in emphysema to macroscopic decline in stiffness using a 3D computational model, PLOS Comput. Biol. 7 (2011) 4. [23] J.B. West, Regional differences in the lung, Chest 74 (4) (1978) 426–437. [24] E. Roan, C.M. Waters, What do we know about mechanical strain in lung alveoli? Am. J. Physiol.-Lung Cell. Mol. Physiol. 301 (5) (2011) L625–L635. [25] B. Suki, D. Stamenovic, R. Hubmayr, Lung parenchymal mechanics, Compr. Physiol. 1 (3) (2011) 1317–1351. [26] D. Haberthur, et al., Visualization and stereological characterization of individual rat lung acini by high-resolution X-ray tomographic microscopy, J. Appl. Physiol. 115 (9) (2013) 1379–1387. [27] T. Sera, et al., Murine pulmonary acinar mechanics during quasi-static inflation using syncrotron refraction-enhanced CT, J. Appl. Physiol. 115 (2013) 219–228. [28] D.M. Vasilescu, et al., Assessment of morphometry of pulmonary acini in mouse lungs by nondestructive imaging using multiscale microcomputed tomography, Proc. Natl. Acad. Sci. USA 109 (42) (2012) 17105–17110. [29] K.S. Burrowes, et al., Towards a virtual lung: multi-scale, multi-physics modelling of the pulmonary system, Philos. Trans. R. Soc. A-Math. Phys. Eng. Sci. 366 (1879) (2008) 3247–3263. [30] H. Kitaoka, S. Tamura, R. Takaki, A three-dimensional model of the human pulmonary acinus, J. Appl. Physiol. 88 (6) (2000) 2260–2268. [31] A. Okabe, et al., Spatial Tessellations: Concepts and Applications of Voronoi Diagrams. Wiley Series in Probability and Statistics, Wiley, England, 2000. [32] D.S. Johnson, et al., Optimization by simulated annealing—an experimental evaluation.1. graph partitioning, Oper. Res. 37 (6) (1989) 865–892. [33] S. Kirkpatrick, C.D. Gelatt, M.P. Vecchi, Optimization by simulated annealing, Science 220 (4598) (1983) 671–680. [34] A. Leeming, R. Schroter, A model morphology of the pulmonary acinus, Proc. Inst. Mech. Eng. Part H-J. Eng. Med. 222 (H4) (2008) 429–437. [35] C.R. Ethier, C.A. Simmons, Introductory Biomechanics From Cell to Oranisms, Cambridge University Press, New Yourk, 2007. [36] C.B. Barber, D.P. Dobkin, H. Huhdanpaa, The Quickhull algorithm for convex hulls, ACM Trans. Math. Softw. 22 (4) (1996) 469–483. [37] L. Xiao, et al., A semiautomatic segmentation algorithm for extracting the complete structure of acini from synchrotron micro-CT images, Comput. Math. Methods Med. 575086 (2013) 1–10. [38] J.E. Hansen, E.P. Ampaya, Human air space shapes, sizes, areas, and volumes, J. Appl. Physiol. 38 (6) (1975) 990–995. [39] J.E. Hansen, et al., Branching pattern of airways and air spaces of a single human terminal bronchiole, J. Appl. Physiol. 38 (6) (1975) 983–989. [40] R.R. Mercer, J.M. Laco, J.D. Crapo, 3-Dimensional reconstruction of alveoli in the rat lung for pressure-volume relationships, J. Appl. Physiol. 62 (4) (1987) 1480–1487. [41] R.R. Mercer, M.L. Russell, J.D. Crapo, Alveolar septal structure in different species, J. Appl. Physiol. 77 (3) (1994) 1060–1066. [42] H. Parameswaran, et al., Three-dimensional measurement of alveolar airspace volumes in normal and emphysematous lungs using micro-CT, J. Appl. Physiol. 107 (2) (2009) 583–592.
©2015 Elsevier
