2.1.

Tissues

All lung tissues were obtained from lung transplant recipients at UW Hospital Madison, Wisconsin, under a current IRB approved protocol. The normal tissues were from pathologist-defined normal adjacent tissue from biopsies of patients without fibrotic lung disease. Tissues were fixed in formalin and sectioned using a vibratome (Leica VT1200) to $\sim 150\mu \mathrm{m}$ thickness. After sectioning, the tissues were stored in phosphate buffered saline (PBS) at 4°C until they were imaged. During imaging, they were mounted on glass slides in PBS with #1.5 coverslips and Vaseline to seal the slide while imaging. A total of six normal and three IPF independent patient samples were prepared and imaged.

2.2.

Microscope System

The imaging system has been described in detail elsewhere³³ and is only described briefly here. The instrument is built around an upright microscope stand (BX61, Olympus, Center Valley, Pennsylvania) with a laser scanning unit (Fluoview 300; Olympus) that is coupled to a mode-locked titanium sapphire femtosecond laser (Mira; Coherent, Santa Barbara, California). All imagings (SHG and TPEF) were performed with an excitation wavelength of 890 nm and an average power of $\sim 20\mathrm{mW}$ at the specimen using a water immersion $40\times 0.8\mathrm{NA}$ objective. This configuration resulted in lateral and axial resolutions of $\sim 0.7$ and 2.5 microns, respectively. Circular polarization at the focus was used to equally probe all fiber orientations. The microscope has two channels with identically calibrated detectors (7421 GaAsP photon counting modules; Hamamatsu, Hamamatsu, Japan). The two channels permitted simultaneous collection of the SHG wavelength (445 nm) in the forward channel using a 20-nm bandpass filter (Semrock, Lake Forest, Illinois) and the elastin autofluorescence signal in the epifluorescence channel using a 22-nm bandpass filter (583 nm; Semrock).

2.3.

Wavelet/PCA/KNN Analysis

We used a wavelet transform to obtain texture features with PCA and KNN analysis for our classification system.^34–36 The wavelet transform decomposition provides both spatial and frequency domain information, which is intricately related to the scale and orientation of the texture features we seek to characterize in the image data. In this process, the wavelet function is placed on a specific location on the image to determine the correlation coefficients between this function and the local morphology. At that location, the shape of the wavelet function is then anisotropically scaled in two dimensions, which then captures (through correlation) both the width and the orientation of the fibers. This process is then translated to different regions in the image and the local wavelet coefficients are calculated. A pictorial diagram of the process going from raw single optical sections to wavelet coefficients is shown in Fig. 1.

Fig. 1

Flow chart of wavelet/principal component analysis (PCA)/K-nearest-neighbor (KNN) algorithm, beginning with the raw image data, calculating wavelet coefficients, and then performing classification using KNN analysis of the extracted PCAs. DWT, discrete wavelet transform.

Since the study size of six IPF and three normal lungs is small, working with the full set of wavelet coefficients (which characterize the input image in terms of the chosen wavelet basis, here the nine filter Daubechies basis) is problematic. In particular, the number of wavelet coefficients one chooses directly corresponds to the dimensionality of the statistical inference problem that needs to be solved in downstream analysis. If the dimensionality of this space is large, one invariably needs to provide the model with a larger number of images to make the inference well-posed. The solution to this problem is to instead analyze the distribution of the wavelet coefficients in terms of their projections on the principal components (PCs). This corresponds to the axes that explain the maximum variance, describing the full set of images with a low-dimensional representation that is more amenable to traditional statistical analysis.

Once these PCs are obtained (via the covariance matrix of the wavelet coefficient distribution), we set up a machine learning task, which constitutes two main steps. First, we use a set of training images, where the class labels of the images are known to learn the pattern that best distinguishes one group from the other (in a space defined by treating the principal axes as the basis). Second, this pattern is used to classify test images whose class label is not known. For classification, we use a simple KNN classifier, a nonparametric method that works under the assumption that the class of each example is similar to the class of its neighbors in the space of PCA axes (see Fig. 1). In other words, for each test image, we consider the majority of votes of its neighbors, which determines the class label of the test image.

The images used in the analysis all came from three IPF and six normal patient samples that were available to us. Each patient sample had several imaging locations, providing different optical stacks. Then 15 individual optical sections were selected from the middle 60% from each optical stack. The middle regions were chosen to avoid any edge effects where the surfaces can be uneven, and also to avoid any effects of attenuation on the signal intensity. As a result of the different numbers of normal and diseased patient samples, there were 270 IPF and 495 normal available optical sections for a total of 765 images. For the wavelet/PCA/KNN analysis, the sample size of the normal and IPF tissues were size-matched, where 270 of the 495 normal images were randomly selected for a total of 540 images (i.e., 270 IPF and 270 normal). PCA dimensions resulting from the wavelet transform of 540 images were randomized and partitioned into 10 subgroups, each with 54 images for KNN classification and cross-validation. Ten KNN cross-validation trials were run, in which nine groups served as the training set and one group was the testing set; each subgroup served as the test group once, as is common in cross-validation experiments. The MATLAB^® code is freely available upon request.

3.1.

SHG Imaging of Normal and IPF Large Airway and Parenchymal Tissues

The ECM structures of normal and IPF lung display significant visual morphology differences in both the large airway spaces and parenchymal areas (Fig. 2). Figures 1(a) and 1(c) are representative SHG optical sections of the collagen architecture surrounding large airways of IPF and normal lungs, respectively. Figures 2(b) and 2(d) are representative optical sections of the collagen architecture of parenchymal areas of IPF and normal lungs, respectively. By visual inspection, the collagen in the diseased lungs in both the large airway and parenchymal regions is packed into denser regions than in the normal tissue. Images at $10\times$ magnification were also acquired and were not visually different than smaller fields of view.

Fig. 2

Panels are representative single second harmonic generation (SHG) optical sections of (a) large airway of an idiopathic fibrosis (IPF) lung, (b) parenchymal region of IPF lung, (c) large airway of normal lung, and (d) parenchymal region of normal lung. Field of view is $180\mu \mathrm{m}$.

For translational purposes, we need to develop objective quantitative methods. The simplest approach is to apply a threshold and calculate the average pixel intensity and collagen area covered, as has been previously reported.²⁶ The threshold used to eliminate the background signal was determined by measuring the background signal in 15 different locations per image stack and finding the average plus the standard deviation. The applied threshold level was specific to each optical section within the image stack and was determined at each optical section to account for signal attenuation at increasing depths.

This approach showed that there were no statistical differences in SHG intensity between either normal and IPF parenchyma or large airway. The area covered was statistically different between normal and IPF parenchyma ($p=0.008$), where the latter was higher, as might be expected for fibrosis. However, there were no differences in coverage between normal and IPF large airway.

3.2.

Wavelet/PCA/KNN Classification of SHG Normal and IPF Tissues

The largely insignificant results in the previous section demonstrate the need for more in-depth quantitative image analysis and classification. The wavelet transform uses the edges of the collagen fibers to provide quantification of the qualitative differences our eyes naturally detect, providing a robust means independent of human visual biases for classifying tissues. The wavelet/PCA/KNN classifier was developed (Sec. 2.3) for this purpose, which reliably differentiates the diseased from normal lung tissues. In this analysis, both large airway and parenchymal regions of the lungs were combined for the classification. Figure 3 shows plots of the receiver operator characteristic (ROC) curves (true positive versus false positive) for the classification of IPF and normal lung tissues for a few combinations of different PCA dimensions and KNN to demonstrate the dependence on PCA and KNN parameters for the classification of tissue. An area under the curve (AUC) of 1.0 is a perfect classification, where 0.5 is a random outcome and provides no discrimination. In practice, values $>0.9$ suggest excellent test accuracy in clinical applications. Table 1 lists the area under the ROC curve for all the combinations of PCA dimensions and KNN applied to classify the images.

Fig. 3

Receiver operator characteristic curves for different combinations of PCA dimensions and KNN used to classify the normal and IPF tissues: (a) varying the number of PCA dimensions from 5 to 300 while constraining $\mathrm{KNN}=5$, (b) varying the number of PCA dimensions from 5 to 300 while constraining $\mathrm{KNN}=25$, and (c) varying the number of KNN from 5 to 15 while constraining PCA dimensions to 100.

Table 1

Area under the receiver operator characteristic curves for combinations of principal component analysis (PCA) dimensions and nearest neighbors.

The optimal classification was obtained using 20 PCA dimensions and 5 NNs, with a resulting AUC of 0.998. The accuracy of classification is similar using five NNs at all PCA dimensions. As the number of NNs is increased, the accuracy of classification still remains high even at low PCA dimension, but the accuracy increases as the PCA dimensions are increased, noting that as more neighbors are included, the slight changes between the higher-order PCA dimensions allow more accurate classification. We found that all combinations of PCA dimensions and KNN provided excellent classification since the worst obtained accuracy was 94%. In general, it is desirable to use as few PCA dimensions and NNs as possible to avoid overfitting errors, and this is enabled here due to the significant change in collagen fibrillar morphology.

3.3.

Determination of Collagen/Elastin Balance in Normal and IPF Tissues

Elastic fiber formation is also increased in IPF (Ref. 10) and the elastin/collagen ratio may be impacted during disease progression (Ref. 9). Initially, in IPF patients, there is an increase in the collagen deposition, but late-stage IPF is described as having an increased presence of elastin. This balance is important in determining the mechanical properties of the lung matrix, such as ECM stiffness and associated elastic recoil forces. We specifically probed both the collagen (SHG) and elastin (TPEF), where these contrasts were simultaneously excited at the same wavelength (890 nm) and spectrally isolated in separate channels. As the elastin contrast is linearly proportional to the concentration, and SHG is a merged effect of the square of the collagen concentration and its organization, it is not possible to determine their actual molecular ratios. The collected signal of both the collagen and elastin signals is further confounded by the different scattering phenomena when imaging relatively thick tissue samples. However, we analyzed the volumetric ratio of elastin and collagen using the well-documented method: $[E_V-C_V]/[E_V+C_V]$, where $E_V$ and $C_V$ represent elastin and collagen voxel volumes, respectively,^31,37 where the limiting cases of this ratio of $+1$ and $-1$ correspond to all elastin or collagen in the pixel, respectively.

Volume fraction estimates were completed on all imaging stacks taken from the parenchyma of three normal and four IPF patient samples. The TPEF spectrum of elastin is broad and overlaps with other autofluorescence signals (e.g., cellular); therefore, segmentation was required to spatially isolate the elastin signal. The intensity of the autofluorescence cellular components was much weaker than that of elastin, which allows thresholding for successful spatial separation as seen in raw [Fig. 4(a)] and segmented [Fig. 4(b)] images for normal, and raw [Fig. 4(c)] and segmented [Fig. 4(d)] images of IPF tissues. The threshold was set using the average signal intensity of the cellular components for each stack, where all the pixels with gray levels above the threshold value were summed to calculate the volumetric fractional coverage of elastin. Similarly, the background signal of the SHG images was eliminated, allowing all the pixels with a signal above the threshold to be summed to calculate the volumetric fractional coverage of collagen.

Fig. 4

(a) and (c) are representative individual two-photon excited fluorescence optical sections of a normal and IPF diseased lung, respectively. The triangular arrows delineate cells, and the arrow indicates mucous, which both are removed via thresholding techniques. The corresponding resulting segmented elastin images are shown in (b) and (d).

Two representative background-corrected color images of the SHG (green) and elastin (blue) TPEF for normal and IPF tissues are shown in Figs. 5(a) and 5(b), respectively, where the organization of the collagen and elastin are both dramatically different in these cases. Specifically, the elastin in the normal tissues is organized within the confines of collagen, whereas in the diseased tissue, the elastin is disorganized and not exclusively intermingled with the collagen fibers. The elastin/collagen index derived from all parenchymal imaging stacks in normal and IPF tissues is shown in Fig. 5(c), where the resulting ratios were $-0.48$ and $-0.63$, respectively, where these values were statistically significant at the $p=0.07$ level using a student’s paired $t$ test. The normal tissues were more elastic relative to collagen than the diseased tissues, indicative of an altered composition of the matrix. This finding, determined by optical microscopy is consistent with known mechanical consequences of the disease. We note that there was no discrimination of large airways through this method.

Fig. 5

Single optical sections of (a) normal and (b) IPF parenchyma, where blue and green is the two-photon excited autofluorescence from elastin and SHG from collagen, respectively. Field of view is $180\mu \mathrm{m}$. (c) Averaged collagen/elastin ratio $[E_V-C_V]/[E_V+C_V]$ of normal ($-0.48$) and IPF ($-0.63$) parenchymal tissues where limiting values of $\pm 1$ are indicative of all elastin and all collagen, respectively. Normal and IPF parenchyma are statistically different ($p=0.07$).