Abstract
Objective. To define the sensitivity of microcomputed tomography- (micro-CT-) derived descriptors for the quantification of lung damage caused by elastase instillation. Materials and Methods. The lungs of 30 elastase treated and 30 control A/J mice were analyzed 1, 6, 12, and 24 hours and 7 and 17 days after elastase instillation using (i) breath-hold-gated micro-CT, (ii) pulmonary function tests (PFTs), (iii) RT-PCR for RNA cytokine expression, and (iv) histomorphometry. For the latter, an automatic, parallel software toolset was implemented that computes the airspace enlargement descriptors: mean linear intercept and weighted means of airspace diameters (, , and ). A Support Vector Classifier was trained and tested based on three nonhistological descriptors using as ground truth. Results. detected statistically significant differences between the groups at all time points. Furthermore, at 1 hour (24 hours) was significantly lower than at 24 hours (7 days). The classifier trained on the micro-CT-derived descriptors achieves an area under the curve (AUC) of 0.95 well above the others (PFTS AUC = 0.71; cytokine AUC = 0.88). Conclusion. Micro-CT-derived descriptors are more sensitive than the other methods compared, to detect in vivo early signs of the disease.
1. Introduction
Chronic obstructive pulmonary disease (COPD) is a complex and heterogeneous disease. Traditionally, two phenotypes of COPD have been described: obstructive bronchitis and pulmonary emphysema. Because of current smoking trends and progressive aging of the world population, an increase in COPD prevalence and related mortality is expected in the coming decades [1]. Emphysema is defined pathologically as the permanent enlargement of the airspaces distal to the terminal bronchioles, accompanied by destruction of their walls, without obvious fibrosis [2, 3]. At the molecular level, emphysema is an inflammatory-driven process caused by the enzymatic destruction of lung elastin and collagen by neutrophil and macrophage elastase [4]. The process is in most cases induced by cigarette smoking [5].
Animal models are key tools to study the disease. Probably the most widely extended one is the mouse model of elastase-induced emphysema, due to its simplicity and low cost [6]. Efficient and sensitive quantification of the injury is needed for the characterization of the disease and the assessment of therapeutic interventions on this model. In this paper, we describe efficient and sensitive methods to quantify histologically airspace enlargement ex vivo. This can in turn be used as an appropriate reference gold standard for the comparison and evaluation of descriptors extracted from micro-CT and other non-radiological in vivo techniques.
The mean linear intercept () [7] is a stereological metric established in the early sixties, commonly used to quantify emphysema in histological samples. is perceived as an index of airspace size, although formally is a measurement of surface area-to-volume ratio. Computing requires counting linear intercepts (see Figure 1), which is a relatively simple manual task [8]. Nevertheless, it cannot be equated with calculating the airspace size without knowledge of the shapes. Weibel et al. [9] illustrate this problem by comparing the calculated both in a sphere and an ellipsoid of equal volumes. Both objects differ in their surface area, which is larger in the case of the ellipsoid. In consequence, measuring linear intercepts from random cross sections of both objects will yield a volume () to surface () ratio for the ellipsoid lower than for the sphere. Indeed, underestimates the severity of emphysema in heterogeneous samples (i.e., samples containing many small airspaces surrounding a few enlarged ones).
To correct this problem, while taking advantage of the computational power of modern computers, Parameswaran et al. [10] presented three alternative measurements (, , and ) based on the moments of the airspace equivalent diameters. When high moments are used, the largest airspaces are weighted more heavily than smaller ones. Therefore, should be more accurate and robust than to quantify heterogeneous, mild emphysema. The authors tested their hypothesis (i.e., shape dependence of and its inability to detect airspace enlargement on heterogeneous samples) on both synthetic images and a few lung parenchyma images. Following this approach, Jacob et al. in [11] showed that finds significantly more pronounced separation between the control and smoke-exposed mice (2–4 cigarettes/day, 6 days/week, 24 weeks) than .
Here, we present a fully automatic and parallel toolset that calculates the above-described metrics (, , , and ) on entire histological sections at reasonable speed. We then used this in-house-developed software to quantify the emphysema in a large study of elastase-induced emphysema mouse model [12], starting at very early stages.
Although histomorphometry is considered the most accurate method for morphological disease characterization, it is an ex vivo technique. There are several other techniques that can be used to characterize lung injury and emphysema in mice in vivo. This is of tremendous value to assess, for instance, the efficacy of new drugs or therapeutic interventions. Among them, micro-computed X-ray tomography (micro-CT) is an especially appropriate modality for lung imaging (Table 1) [13]. Using micro-CT, moderate to severe forms of emphysema can be easily detected in elastase-treated mice, as the lung parenchyma appears darker than in nontreated mice (see Figure 2). Pulmonary function tests (PFTs) provide physiological data about lung function and their characteristic parameters (Resistance and Compliance) [12]. Finally, the inflammatory response of the lungs to the damage caused by the insult can be evaluated by looking at protein or RNA cytokine levels, using blood plasma or lung tissue extracts. The information provided by these methods is complementary.
(a)
(b)
In this study, we train a classifier to evaluate and compare the sensitivity of micro-CT density-based-derived descriptors and other non-radiological techniques (i.e., single compartment model PFTs, and RNA inflammatory cytokine levels) to detect elastase-induced lung damage. The histomorphometrical parameter is used as reference value because is a generalized, shape-independent quantitative descriptor of airspace dimensions, well suited for automated computation.
The paper structure is as follows. In the Materials and Methods, we describe the animal preparation, the PFTs and micro-CT imaging, and the sample preparation for RNA cytokine levels measurement and histology. We also present the image acquisition, the image analysis pipeline, and the statistical analysis of the histology data. Besides, we describe the classifiers that we used and how they were trained. The next section presents the results on the sensitivity of the histomorphometric descriptors to detect differences between the control and elastase-treated groups a few hours after the insult, and to study disease progression. Finally, we report on the ability of the rest of descriptors to adequately classify the animals using the histology as reference. The paper ends with a discussion and conclusions.
2. Materials and Methods
2.1. Animal Preparation and In Vivo Tests
2.1.1. Animal Preparation
All experimental protocols involving animal manipulation were approved by the University of Navarra Experimentation Ethics Committee. Sixty, 11-week-old mice were equally distributed into a control and a treatment groups. Treated mice were intratracheally instilled with 6 units per 30 g of porcine pancreatic elastase (PPE, EC134GI, EPC, MI, USA), as described in a previously published protocol [14]. Control animals were instilled with a saline solution. Five animals of each group were sampled at hours 1, 6, 12, and 24 and days 7 and 17 after treatment. At each time point, all animals underwent micro-CT thoracic imaging and pulmonary function tests. Then, they were sacrificed and their lungs collected for histomorphometry and cytokine measurements (RNA and protein).
2.1.2. Pulmonary Function Tests
Animals were anesthetized, intratracheally cannulated, and connected to a Flexivent rodent ventilator set (Scireq, Montreal, QC, Canada) set at a rate of 200 breaths/min and a tidal volume of 10 mL/kg. Lung resistance () and compliance () values were measured using a single-frequency-forced oscillation, fitting the measured data to a single compartment model [15]. All measurements were repeated three times.
2.1.3. Breath-Hold-Gated Micro-CT Imaging
The anesthetized, artificially ventilated mice were scanned using an X-ray micro-computed tomograph (Micro-CAT II, Siemens Pre-Clinical Solutions, Knoxville, TN, USA). Parameters used for the micro-CT image acquisition are those given in Table 2. Seven hundred micro-CT projections were acquired during 650 ms isopressure breath holds at 12 cm H2O. Two normal breathing cycles were induced between breath holds. A total lung capacity perturbation was performed every 20 breath holds to prevent atelectasis.
The reconstructed three-dimensional images had a total of 640 slices with isotropic 46 μm voxel size and pixels per slice. The scanning time was approximately 30 minutes. The estimated X-ray dosage was 71.6 cGy/exam (Siemens Pre-Clinical Solutions, Knoxville, TN, USA). A water phantom was used to calibrate the image to Hounsfield units (HU).
2.1.4. Micro-CT Image Analysis
First, the lungs were automatically delineated in the 3D micro-CT images using an existing segmentation method [16]. Then, the airways were segmented and reconstructed using a fast and robust algorithm developed by us [17]. The algorithm is based on a propagating fast marching wavefront that divides the tree into segments as it grows. A number of specific rules were applied after every iteration of the front propagation to avoid unwanted leakage into the parenchyma. From the segmented airways, the right and left radius measurements of the mainstem bronchi (RMBR and LMBR, resp.) were computed. The airways were removed from the lung volume before quantification. Two other parenchymal injury descriptors were also computed: mean lung voxel intensity (MLVI) and relative volume below −900 HU (VBT). This threshold was selected because intensity values below −900 HU are rare in scans of healthy mice (the VBT is less than 5% of the total lung volume in all healthy animals of any age).
2.2. Sample Preparation for RNA Cytokine Expression and Histology
Mice were anesthetized and then sacrificed by exsanguination. Next, the lungs were fixed at a constant pressure of 20 cm H2O. RNA was extracted from the accessory lobe and then qRT-PCR was performed with an Applied Biosystems 7900HT Fast Real-Time PCR System. Four immune-modulatory cytokines and chemokines involved in lung inflammatory response [18] were analyzed, namely, interleukin 6 (IL6), immune protein 10 (IP10), keratinocyte chemoattractant (KC), and monocyte chemoattractant protein 1 (MCP1). Beta-2 microglobulin (B2M) was used as the endogenous control gene for the experiments.
Three paraffin blocks from different lobules of each lung were created and reserved for histological analysis. Three nonconsecutive sections per block were cut and stained with Hematoxylin and eosin (H&E), resulting in a total of nine slides per mouse, each containing two lobe pieces. The total number of sections acquired and analyzed was 1080.
2.3. Histological Image Analysis
2.3.1. Image Acquisition
Whole-slide views of all 1080 sections were acquired using an automated Axioplan 2ie Zeiss microscope (Carl Zeiss, Jena, Germany). Each slide was initially acquired with a Plan-Neofluar objective (numerical aperture NA = 0.035, magnification 1.25x, pixel resolution 3.546 μm/pixel). The automatic threshold method proposed by Otsu [19] was then applied to detect all tissue areas. The size of the objects was measured and only objects with a reasonable size to represent entire sections of lung lobes were considered for further processing. For each object, a bounding box was created and the coordinates of its four vertices were sent to the microscope. Then, an automatic routine scanned those areas with a Plan-Neofluar objective (NA = 0.3, 10x, 0.725 μm/pixel). Some overlap was allowed between image fields to facilitate the creation of large mosaics. The Stitcher ImageJ plugin [20] was used for it. The resulting mosaics were stored in a server for quantitative analysis. Figure 3 shows a sample mosaic (Figure 3(a)) with a zoomed-in area showing heterogeneously distributed airspace enlargement (Figure 3(b)).
(a)
(b)
2.3.2. Image Analysis
To accurately quantify bona-fide emphysematous air spaces, all vessels, alveoli, and spurious structures must be recognized and removed from the images because they might be confused with enlarged alveolar spaces. In the following paragraphs, we describe the method used to segment vessels and alveoli and the extraction of image descriptors.
Segmentation
The main steps of the segmentation algorithm are summarized in the flowchart shown in Figure 4 and illustrated the snapshot shown in Figure 5. This image contains a vessel of a certain size in the center of the field of view. First, the 8-bit grayscale green channel was extracted from the 24-bit RGB, since it provides the greatest contrast between the background and the red-blue H&E stained tissue. Then, a mask of the parenchyma tissue was obtained by thresholding the histogram of the image. The histogram of a typical histological image is monomodal. Thus, common bimodal thresholding techniques—such as Otsu’s thresholding—would not work well, leaving out some tissue structure walls. As an alternative, a maximum deviation unimodal thresholding [21] was used (Figure 6(a)) to extract all tissue areas from the background. The obtained mask was then inverted and all the luminal structures were segmented and labeled by connected component labeling (Figure 6(b)).
Since the lumen of blood vessels and bronchioles could be confused with enlarged airspaces, they must be detected and removed from the image before proceeding with the quantification. To this end, we applied binary erosion to the mask of the parenchyma using a structuring element of size 7 to remove all but the thickest walls, corresponding to vessels and bronchioles. Then, we used the geometric methods presented by [22, 23] and implemented in the Shapely Python package [24] to calculate the convex hull of the remaining walls (see Figure 6(c)). Next, the intersection between the convex hulls of the walls and the centroids of the labeled structures is performed and all the structures for which a nonempty intersection exists (i.e., vessels and bronchioles) are removed from the segmentation (see Figure 6(d)). At the end of this process, all the labeled regions represent airspaces.
The computation of the intersection is an intensive operation. To increase the speed, this operation was performed on a downsampled version of the image. The downsampling factor was empirically set to four to avoid loosing wall structures. However, the extraction of the airspace enlargement descriptors is performed on the full-resolution image.
(a)
(b)
(c)
(d)
Extraction of Descriptors
We extract four descriptors from the histology, namely, the linear mean intercept and the moments of the airspace equivalent diameters (, , and ).
The classical mean linear intercept is defined as the mean length of the linear intercepts in the lung. In this paper, an approximation of is computed as follows: the labeled image (Figure 6(d)) of each lobe is converted into an array of pixels; the image is then raster scanned, and the number of positive pixels between each pair of two consecutive null (wall) pixels is counted and its value (i.e., length of the th linear intercept ) is stored; then, the is calculated as the mean of all the stored linear intercepts lengths on the lung lobe.
To compute the indexes, we first calculate the area of the th airspaces () of the segmented, labeled image, by counting the number of pixels inside each labeled region. This value is scaled to physical dimensions using the objective calibration.
Then the equivalent airspace diameter of each space is defined as
which is equal to the diameter of a circle with area .
The family of indexes is then defined as the ratio of two moments of the equivalent airspace diameters distribution :
where indicates the arithmetic mean. can be expressed as functions of the central moments of airspace diameters. In particular, is the arithmetic mean of the airspace diameters. is a function of its mean () and variance ():
and is a function of , , and skewness () of the airspace diameters:
Computational Architecture
The code was written on a hybrid Python/C++ platform and was specifically designed to exploit parallelism. First, a process is spawn to walk the image directory tree and collect the files to be processed. Then, a job description is created for each image, which contains all the parameters used to guide the analysis through the steps detailed before. Each job is added to a shared queue on the network. Several processes access the remote queue from five different machines, download a job description, and execute it locally. Results are saved in a distributed file system.
2.3.3. Statistics
The median and interquartile ranges (IQR) were calculated for each histology descriptor and time point. The control and elastase-treated groups were compared using the Mann-Whitney U-test. To measure progression, the same test was performed on the elastase-treated animals at successive time points. values ≤ 0.01 were considered statistically significant. The language and environment for statistical computing [25] were used for statistical analysis.
2.4. Classification
In the previous subsections, we have presented three different sets of disease descriptors: micro-CT density-based descriptors: MLI (HU), VBT (%), RMBR (μm), and LMBR (μm); Pulmonary Functional Test single compartment model parameters: Resistance () (cm H2O/mL/s) and Compliance () (mL/cm H2O); RNA cytokine expression measured as concentration: IL6, IP10, KC, and MCP1.
A Support Vector Machine (SVM) Classifier was created for each descriptor set. The chosen kernel was a Gaussian radial basis function. The classifier was defined by two parameters: the area of influence of the support vector on the data space and the soft margin parameter . The best parameter combination was selected by a grid search with exponentially growing sequences of and . In particular, was chosen in the interval (10−9, 103) and in (10−1, 1013). Each combination of parameters was cross-validated, and the parameters with best cross-validation accuracy were selected. The dataset was divided in a training (60% of the samples) and test groups (the remaining 40%).
For every classifier, a receiver operating characteristic (ROC) curve was generated. The histomorphometrical parameter was used as the gold standard and the 3rd quartile value of the control animals was used as a threshold. For all possible combination of features, the classical performance indexes of Area Under the ROC Curve (AUC) and -score on the test dataset were computed. The results were then compared to determine the best classifier.
3. Results
3.1. Histological Analysis
The mean linear intercept () and weighted mean (, , and ) were computed and used as estimates of the airspace enlargement. Figure 7 and Tables 2 (), 3 (), 4 (), and 5 () show the progression of these parameters throughout the experiment for both groups (control and elastase induced). As expected, no progression was found in the control group, regardless of the parameter used. In the treatment group, all the parameters increased significantly () from 1 hour to 24 hours after elastase treatment. and have also a statistically significant increase () from 24 hours to one week after elastase treatment, which was somehow undetected by and . No progression was detected in any case from one week to 17 days after the insult.
(a)
(b)
(c)
(d)
Regarding the sensibility to differences between the control and treatment groups, is significantly different at all time points except at 6 and 12 hours after elastase administration. is significantly different at all time points except very early on, 1 hour after treatment. Contrarily, and are significantly different () between treatment and control groups at all time points. In summary, although all the parameters yield similar results, and are more sensitive to differences between the emphysematous mice and controls.
The computation time per lung lobe section was about 20 minutes.
3.2. Classification
We used three classifiers based, respectively, on: Micro-CT density-based descriptors, RNA cytokine expression measured as relative RNA concentration by qPCR, and single compartment model parameters from pulmonary function tests.
The cross-validation accuracy plots for the selection of the optimal estimated parameters are shown in Figures 8, 9(a) and 9(b) for the micro-CT, cytokine expression and pulmonary function tests classifier, respectively. The optimal estimated parameters, best set of features, AUCs, and -scores are given in Table 6 while Figure 10 shows the ROC curves for the best classifier of each parameter set. The Micro-CT-based classifier using as features MLI, VBT, and RMBR shows the best performance with an AUC of 0.95 and a -score of 0.92 and as shown by the ROC curve the highest true positive rates can be achieved for the same false positive rate. In general, RNA cytokine expression classifier using as features KC, IL6, and IP10 performs slightly worse than micro-CT with an AUC of 0.88 and -score of 0.71 although it achieves the best true positive rate of 0.8 for a false positive rate of 0.1. The classifier based on feature Functional resistance performs significantly worse than the other two classifiers with an AUC of 0.71 and an -score of 0.66, being the closest to the performance of a random classifier.
(a)
(b)
4. Discussion
The main aim of this work was to evaluate and compare the sensitivity of the quantification of elastase-induced lung damage, using micro-CT-derived descriptors, pulmonary function tests based on a single compartment model, and RNA cytokine expression. Histomorphometry was used as gold standard for the comparison.
Our results show that is able to distinguish between the control and the elastase-treated group at all time points, starting as early as one hour after treatment. This early airspace enlargement might be caused by surfactant dysfunction resulting from elastase administration. The ability of to discriminate such early damage can be attributed to the fact that heavily weights on enlarged airspaces and therefore, reflects better the airspace size distribution. In terms of the disease progression, and detected airspace enlargement during the first 24 hours and from that point until one week after treatment. This last increase was missed by and .
Compared to previous studies, our histomorphometry values were obtained using a considerably larger sample size. In previous works, a few random fields were acquired from each mice lung. Here instead, mosaic images of whole lung lobe sections were acquired and analyzed, which was possible thanks to our fully automated software. Other interesting characteristics of automated systems like ours are the reduction of operator bias and a substantial reduction of the time dedicated to the analysis. Finally, the key for the success of our method is the fact that our toolset implements a mechanism to eliminate alveoli and vessels. Previously, major airways and vasculature were manually discarded at acquisition or analysis time using a great deal of manual interaction.
Micro-CT is especially appropriate to study in vivo the progression of animal models of pulmonary disease [13]. We have set up a generic protocol for micro-CT image acquisition that allows longitudinal studies [26]. The protocol includes endotracheal intubation and iso-pressure breath holds to reduce movement artifacts. Several segmentation and analysis methods were developed to quantify the effects of disease on the very noisy, artifact-plagued micro-CT images [12, 16]. These methods allow for quantitative measurements of the lungs and the airways separately, thus allowing to monitor disease development. In this work, we found that a SVM classifier using the micro-CT-derived, features MLI, VBT, and RMBR reached a high AUC and -score, thus indicating that micro-CT produces reliable measurements of airspace enlargement even at very early disease stages. Pulmonary function tests were also performed using forced oscillation techniques with endotracheal intubation. The SVM classifier trained on the Pulmonary function tests parameters using only tissue resistance () achieved the best AUC and -score. It could seem strange that the optimal classifier uses instead of or both. In long-term studies of the elastase-induced model, gets clearly increased as a reflection of high lung stiffness. Our results presented in [12] show instead that decreases during the first 24 hours. It is only at week 4 that starts to increase. We hypothesize that this behavior may be related to the acute inflammatory reaction occurring in the elastase-induced model, starting immediately after treatment and nearly disappearing by day 7. On the other hand, the relatively good performance of the SVM classifier training on the cytokine expression levels could reflect this short-time inflammation. However, its ability to detect airspace enlargement in long-term studies has yet to be confirmed.
Finally, as explained in detail in the Results, the micro-CT-derived descriptors seem to be especially well suited for the estimation of airspace enlargement on the elastase-induced emphysema mice model.
5. Conclusion
In this paper, we presented open source software for the automatic quantification of airspace enlargement in large histological tissue sections. Using this software, we can automatically process large amounts of data in a relatively short period of time and with minimal user interaction. The automated measurements were able to detect airspace enlargement very early after the insult with elastase. Those measurements were used as ground truth to assess the sensitivity of micro-CT and other non-radiological techniques. Interestingly, typical respiratory-gated micro-CT density-based descriptors (mean lung density and relative volume below −900 HU) and the right radius of the mainstem bronchi achieved a high sensitivity and specificity discriminating early disease signs.
Abbreviations
| AUC: | Area under the ROC curve |
| B2M: | Beta-2 microglobulin |
| : | Compliance |
| COPD: | Chronic obstructive pulmonary disease |
| : | Mean alveoli area |
| : | corrected by the first order statistical moment |
| : | corrected by the second order statistical moment |
| H&E: | Hematoxylin and eosin |
| IL6: | Interleukin 6 |
| IP10: | Immune protein 10 |
| KC: | Keratinocyte chemoattractant |
| LMBR: | Left mainstem bronchus |
| : | Linear intercept |
| Micro-CT: | Microcomputed X-ray tomography |
| MLVI: | Mean lung volume intensity |
| MCP1: | Monocyte chemoattractant protein 1 |
| PFT: | Pulmonary function test |
| : | Lung resistance |
| RMBR: | Right mainstem bronchus |
| ROC: | Receiver operating characteristic |
| SVM: | Support vector machine |
| VBT: | Volume below threshold. |
Acknowledgments
This project was partially funded by the “UTE Project CIMA;” the Spanish Ministry of Health under Grants PI070751 and RTICC RD06/0020/0066; the subprogram for unique strategic projects of the Spanish Science and Innovation Ministry under Grants MICINN PSE SINBAD and PSS 0100000-2008-2; the program for nonfundamentally directed research Grants of the Spanish Science and Innovation Ministry under Grants MCYT TEC2005-04732, MICINN DPI2009-14115-C03-03, and MINECO DPI2012-38090-C03-02. A. Muñoz-Barrutia holds a Ramón y Cajal Fellowship of the MICINN. M. Ceresa holds a Torres y Quevedo Fellowship of the MICINN.