2.1.

Animal Models and Experimental Protocols

The animal protocol has been approved by the committee on animal care and use at the University of Maryland, College Park. A murine model of CKD was induced by injection of ADR ($1.5\mathrm{mg}/\mathrm{kg}$) into the tail vein of Munich-Wistar rats.³³ Once a week during the entire study protocol of 8 weeks, rats were weighed. The 24-h urine volumes were collected in metabolic cages. Fresh urine samples were tested for albuminuria (Albustix), and blood samples taken from the tail vein were analyzed for serum creatinine and BUN values (Beckman Coulter Creatinine and BUN Analyzers). In each week for the 8-week period, two rats were anesthetized with $\text{isoflurane}/{\mathrm{O}}_2$ (4% induction, 1.5% during operation, ${\mathrm{O}}_2$ $1\mathrm{L}/\min$). The abdominal cavity was opened through a midline incision, and the left kidney was exposed and imaged using OCT. Following in vivo OCT evaluation, the kidneys were fixed in situ by flushing with warm (i.e., 37°C) oxygenated saline, followed immediately by phosphate buffered 2% paraformaldehyde and 0.1% glutaraldehyde. The fixed kidneys were excised and the rat euthanized by intracardiac injection of pentobarbital sodium to induce cardiac arrest. Blocks of fixed kidneys were embedded in paraffin, sectioned, and stained with hematoxylin and eosin.

2.2.

Optical Coherence Tomography

A custom-built Fourier-domain OCT system was used in this study (Fig. 1).^34–37 This system uses a swept-source laser operating at 1310-nm center wavelength with 100-nm bandwidth. The axial and lateral resolutions of the system are 12 and $6\mu \mathrm{m}$, respectively. OCT image dimensions are 1024 pixels ($X=1.4\mathrm{mm}$) in lateral direction and 512 pixels ($Z=2.0\mathrm{mm}$) in axial direction. The sensitivity of the system is 90 dB. The A-scan acquisition speed is 16 kHz. The 2-D image acquisition speed is 16 frames per second. For each experiment, 3-D OCT volumes were acquired from 5 to 10 kidney locations. Each volume consists of 475 consecutive 2-D images.

Fig. 1

The schematic diagram of the custom-built OCT system used in this study. Inset image shows the abdominal cavity of a Munich-Wistar rat opened through a midline incision and the exposed left kidney for OCT imaging.

2.3.

Computer-Aided Diagnosis Software

About 3000 to 5000 OCT kidney images have been acquired each week for eight weeks. To determine CKD condition from the massive image dataset, an OCT image analysis software is developed to automatically detect and measure features that are related to CKD progression. In our application, possible features can be the proximal tubule morphology, such as tubular diameter and atrophy/hypertrophy.³⁸ The image analysis software consisted of four sections: image preprocessing, feature region proposing, feature region classification, and feature region measurement.

2.4.

Image Preprocessing

Gaussian blur is used to reduce the speckle noise of the OCT image. Gaussian blur is performed by convoluting the original OCT image with a 2-D Gaussian kernel. As a low-pass filter, Gaussian kernel removes high-frequency noise as well as useful high-frequency information, such as edges. A small Gaussian kernel ($\text{sigma}=2$) is used to limit the blur effect and retain small and fast-varying features as much as possible.

2.5.

Feature Region Proposing

To recognize multiple objects in the image, regions of interest (ROIs) have to be proposed to an image classifier. One of the easiest ways to propose ROIs is the sliding window algorithm. Specifying the range of window size and the stride size, a window slides from the edge of the image and proposes each window patch as an ROI. This method has been used widely and successfully in applications such as face recognition.³⁹ However, as a greedy search method, it proposes a large number of proposed regions and suffers from “curse of dimensionality.”⁴⁰ Sliding window method has to be combined with a very fast image classifier for speed-critical applications.

Another region proposing method is segmentation. Segmentation divides the image into subimages based on low-level image properties, such as pixel intensity and texture. For OCT images, image intensity-based segmentation with a global threshold cannot recall tubule lumens with high accuracy due to the nonuniform image brightness. A way to overcome this issue is to remove the effect of illumination in preprocessing, such as by the Retinex theory.⁴¹ However, in practice, those enhancement algorithms are very slow. In addition, the illumination effect may not be removed perfectly. Another way is to use dynamic local intensity thresholding. “Bradley” adaptive thresholding technique uses a computationally efficient integral image algorithm to determine the local threshold value for each subimage window.⁴² In comparison to other dynamic thresholding algorithms, Bradley adaptive thresholding is robust and less computationally expensive.⁴³ A representative result of an OCT image segmented by Bradley adaptive thresholding is shown in Fig. 2(a).

Fig. 2

Same OCT image processed by (a) dynamic local intensity thresholding and (b) SLIC superpixel.

Superpixel is another region proposing method that groups similar pixels into a “superpixel.” Those superpixels can then be proposed as ROIs. The state-of-the-art superpixel methods include graph-based algorithm,⁴⁴ gradient-ascent-based algorithm,⁴⁵ and SLIC algorithm.⁴⁶ A representative result of an OCT image processed by SLIC superpixel method is shown in Fig. 2(b). The principle of SLIC superpixel algorithm is the following: (1) first, $N$ equally spaced pixels are initialized as cluster centers, (2) each pixel then associates itself with the most similar neighboring cluster centers, (3) the cluster centers are updated to be the average pixel location of all its associated pixels, (4) repeat steps 1 to 3 until the new cluster center and the old cluster center converges. The number of cluster centers $N$ controls the granularity of superpixels or the total number of proposed regions.

Ideally, the region proposing algorithm should have a 100% recall rate. Recall rate is defined as the percentage of total tubule lumens being proposed. Another key parameter is the total number of proposed regions. For example, sliding window algorithm with a small stride size and a large range of scale would have a near 100% recall rate, because it can theoretically capture all tubule lumens in different scales and at different locations. However, its total number of proposed regions is massive, which is impractical for many applications.

We compared the ROI recall rate and the total number of proposed regions between Bradley segmentation and SLIC superpixel. To determine the recall rate, tubule cross-sectional regions in test images are manually labeled [Fig. 3(a)]. A tubule region is correctly recalled if more than 60% of the recall region overlaps with the manually labeled region [Fig. 3(b)]. The number of proposed regions can be altered by changing the threshold level and window size for Bradley adaptive threshold algorithm or the initial number of cluster centers for SLIC superpixel algorithm. The recall rate versus the number of proposed regions for each algorithm was measured [Fig. 3(c)]. From the results, Bradley adaptive thresholding with a window size of 30 pixels is able to recall more than 80% of the tubule regions with under 500 proposed regions. SLIC superpixel algorithm underperforms Bradley adaptive thresholding for this application. Therefore, we choose Bradley adaptive thresholding with a window size of 30 pixels in this study.

Fig. 3

(a) Manually labeled tubule lumens (green regions), (b) overlaps between proposed regions (white) and manually labeled true regions (green), and (c) recall rate versus number of proposed regions for different algorithms.

After candidate regions for proximal convolutional tubule lumens are chosen, rectangular image patches containing the candidate regions are cropped for further image classification. It is unavoidable that the set of proposed images would include many false positive regions, such as defects in the OCT image. From the perspective of adaptive threshold, they are similar to tubule cross sections, which have a relative darker center than the surrounding regions. Examples of the region-proposed images are shown in Fig. 4. To differentiate them, an image classifier is necessary to distinguish the true tubule cross sections [Fig. 4(a)] from the falsely proposed regions [Fig. 4(b)].

Fig. 4

Examples of proposed region images: (a) tubule lumens and (b) OCT image artifacts and speckle noises. The brightness of some images has been increased for the ease of viewing. Sizes of the proposed region images range from 5 to 60 pixels in width.

2.6.

Feature Region Classification: Convolutional Neural Network Image Classifier

There are many image classifier models. Some of the well-known classifiers include hard-coded classifier, eigen features classifier,⁴⁷ and artificial neural network,⁴⁸ etc. A classifier generated through learning tends to be more accurate due to the amount of prior information gained through the training process. In our study, OCT images contain large variations in tubule morphology, kidney condition, image quality, etc. ConvNet⁴⁹ is eventually chosen as the image classifier due to its larger learning capacity.

ConvNet training requires thousands to millions of labeled images per detection class depending on the complexity of the network structure. In general, a large and complex network is capable of learning to classify more complicated visual features. However, it would also require a large annotated training database. Overfitting happens if a large network is trained with a small database, which leads to low-classification accuracy in testing. Therefore, there is a trade-off between learning capacity (network size) and training database size. In practice, despite the fact that sufficient image data are usually available, for a supervised learning model, such as ConvNet, it is usually the labor cost for annotating the thousands to millions of images that make creating a large training database a prohibitive task.

For our application, we built a training database with two classes: tubule class and nontubule class. To crop and collect training images, we first created a less-accurate image classifier based on hard-coded features. Bootstrapping from this classifier, we were able to collect sufficient images per class easily but with significant false classification rate. The incorrectly classified images in both classes were manually removed from the database by inspection. To keep the manual inspection time minimal, we deliberately collected only 2000 training images per each class. However, 2000 images per class are far from enough to train a network that has thousands of weights and biases without overfitting. To increase the database size, each of the original images was duplicated $N$ times. Then, the duplicated images went through random translation, rotation, scaling, and noise addition (Fig. 5). The database size was artificially increased by $N$ times. In our case, $N$ equals to 10 and the final database size is 20,000 images per class.

Fig. 5

Multiple training images generated by manipulating the original tubule cross-section image (first on the left). All images are resized to $32\times 32\text{pixels}$ for ConvNets training.

We attempted to train two different ConvNets with the database. The first network has a small structure size (less than the overfitting limit of the database) as shown in Fig. 6(a). The second network is the LeNet⁴⁸ [Fig. 6(b)]. LeNet is chosen, because it is originally designed for recognizing single channel grayscale handwritten-digit images (MNIST dataset), which is in the similar format as OCT images. Many latest ConvNet architectures, such as AlexNet, are originally designed for RGB image classification; therefore, they are not included in this study. The smaller-sized network was trained end to end from the first to the last layer with each layer initialized randomly. LeNet was trained with transfer learning; instead of end-to-end training, it started from a pretrained LeNet for handwritten-digit recognition. Only the fully connected layers were fine-tuned using our training database. Transfer learning enables a large network to be trained on a smaller dataset without overfitting. The performance of the two networks was benchmarked with a separate testing database; the testing database contains OCT images randomly chosen from all CKD stages. Locations of the tubules in the testing database were manually labeled. The bench test result showed a better tubule recognition rate (90%) from LeNET with transfer learning. It is worth noting that whether transfer learning is a good option can be speculated based on the degree of visual differences between the original application of the pretrained network and the current application.⁵⁰ The ConvNet training was performed using MatConvNet⁵¹ on MATLAB^® using an Nvidia GTX660 GPU. Once the network was trained, the network weights were saved. Image classification based on the weights could be performed on machines without Nvidia GPU support. In our study, the image classification based on the trained ConvNet was programmed in C++, with multithreading parallelism using OpenMP, and ran on an Intel 8-core Xeon Processor. The processing speed is about 10 frames per second.

Fig. 6

(a) An end-to-end trained small ConvNet and (b) LeNet with only the fully connected layer being trained while the rest (gray) of the layers inherited from the pretrained network trained on MNIST database.

2.7.

Feature Region Measurement

CKD is defined as the progressive loss of kidney filtration function in time. Nephrons are the basic filtration unit in the kidney. A histology study of a rat CKD model has shown that when nephrons lose their filtration capability, the proximal convolutional tubules of those nephrons shrink and close.⁵² The rest of the functional nephrons in the kidney would compensate the dysfunctional ones. The proximal convolutional tubules of the functional nephrons enlarge and become hypertrophic. A histology study has shown that the average convolutional tubular diameter increases as CKD condition progresses.³³

For each tubule’s cross-sectional region, the circularity of the region was first measured. Since a 2-D cross-sectional OCT image could cut the proximal tubule in any random plane, the circularity of the region indicates whether the tubule cross section is a perpendicular cut. Only the diameters of tubule regions with high circularity ($>0.8$) are chosen for diameter analysis. The procedure to measure the circularity and diameter of a lumen is the following. (1) Locate the contour pixels of the region. The circumference of the region ($L$) is the total number of contour pixels. (2) Locate all the pixels that belong to the region. The area of the region ($A$) is the total number of pixels in the region. (3) Circularity of the region is calculated by dividing the circumference square of the region by $4\pi$ times of region area, i.e., $L^2/(4\pi A)$. A perfect circle would have a circularity equal to 1. (4) Locate the region center by averaging the coordinates of each pixel inside the region. (5) Distances from the region center to each contour pixel are measured. The diameter of the region is 2 times the mean of the distances. Figure 7 shows the process to determine the circularity and diameter of tubule cross-sectional region.

Fig. 7

Illustration of definitions used in region measurements.

3.1.

Optical Coherence Tomography Imaging of Chronic Kidney Disease Rats

Visual differences between normal kidneys and CKD kidneys can be perceived from OCT images of rat kidneys. Images in Fig. 8 are collected from both healthy rats and rats with mid- to late-stage CKD conditions. In healthy kidneys, the proximal tubule lumens are similar in size and distributed uniformly with a homogeneous distribution pattern [Figs. 8(a)–8(c)]. In kidneys with late-stage ADR-induced CKD [Figs. 8(g)–8(i)], some proximal tubules either appear reduced in size or become indiscernible to OCT, i.e., atrophic. At the same time, the rest of the proximal tubules appear increased in size, i.e., hypertrophic. This is consistent with the published histology study in dysfunctional nephrons.³³ These are the two most extreme cases in CKD kidney morphology. In the midstage CKD (week 2 to 6), the observed kidney images [Figs. 8(d)–8(f)] can be a mixture of these two extreme cases, where the hypertrophic tubules are smaller than those at late-stage and atrophic tubules are more discernible than those at late stage.

Fig. 8

(a)–(c) are images captured from healthy kidney, where proximal tubules are uniform in size, (d)–(f) are images captured from midstage CKD kidney, where certain percentages of proximal tubules appear shrunk and certain percentages of proximal tubules appear enlarged, (g)–(i) are images captured from late-stage CKD kidney, where proximal tubules in some cross sections appear further enlarged and in some other cross sections disappeared due to shrunk to extreme small sizes.

3.2.

Individual Image Inspection

While it is possible to determine whether CKD exists, it is difficult to quantify the severity of the disease. To quantify the progression, we developed a CAD software to analyze the OCT images and measure the size of the proximal tubules. The CAD software was first tested on individual OCT images captured from different animals at different stages of CKD (Fig. 9). The software can detect and measure the tubular morphology and present the statistics instantly as the user navigates through images. In addition, the software highlights the detection using color-coded bounding boxes, where yellow, blue, and green indicate normal-size tubules, hypertrophic tubules ($\text{diameter}>50\mu \mathrm{m}$), and atrophic tubules ($\text{diameter}<15\mu \mathrm{m}$), respectively. Figure 9 shows two images examined by the software. Figures 9(a) and 9(b) are from a healthy kidney and a late-stage CKD (8 weeks after ADR injection) kidney, respectively. Visually, it is apparent that tubule lumens in Fig. 9(b) are larger than those from Fig. 9(a). The software measures the tubule lumen diameters, and the results from the software show that tubules from healthy kidneys are mostly within the normal ranges, which is 20 to $30\mu \mathrm{m}$. While the tubules from late-stage CKD kidneys are either in the hypertrophic range ($\text{diameter}>50\mu \mathrm{m}$) or in the high end of normal ranges. The result is consistent with previous reports that the diameters of the hypertrophy tubules can be 50% larger in mice and nearly 100% larger in human than the diameters of the normal tubules.^33,38 These results verify that the CAD software is able to detect and measure tubule cross sections accurately, which can be used to differentiate healthy and nonhealthy kidneys.

Fig. 9

Individual image analysis shows significant tubular morphological differences between (a) healthy kidneys and (b) late-stage CKD kidneys.

3.3.

Automatic Batch Analysis

In addition to individual image inspection, the software is able to analyze the entire batch of image dataset automatically. The processing speed is about 10 images per second. It takes about 10 min to process, measure, and generate statistics of a dataset containing 5000 images. The software records the information of each detected tubule, including diameter and circularity. The OCT image dataset can be quickly analyzed by the software without supervision of the user.

The complete OCT image dataset of rat kidneys was analyzed by the software. The dataset consists of 18 animals in total: 2 animals per week group for 8 post-ADR-injection weeks in total and 2 healthy animals for control. 3000 to 5000 2-D OCT images were collected per rat. In total, the complete dataset contains more than 70,000 images. The CAD software measures and records tubule diameters and generates the statistics, including tubular diameter, percentage of hypertrophic tubules, tubule density, and tubule diameter range for each week.

Tubule diameters are measured directly by the software after the tubule being detected. Percentage of hypertrophic tubules is derived by recording the number of detected tubules whose diameters are larger than $50\mu \mathrm{m}$. Tubule diameter range is the difference of the tubule diameter at 95 percentiles and that at 5 percentiles. Tubule density is computed by determining number of tubules per unit tissue area.

Figure 10 shows representative OCT images for the control and CKD-induced rat. For the healthy animal, most of the tubule diameters are within normal range (yellow). One week after injection, fewer tubules are visible and more atrophic tubules are present. During the midstages (week 2 to 6), number of both atrophic tubules and hypertrophic tubules increases. This is probably due to the fact that functional nephrons overwork to compensate the loss of atrophic tubules (green). At the late stages (week 7 to 8), more hypertrophic tubules (blue) are visible.

Fig. 10

Sample kidney OCT images from healthy control rats and from CKD rats at each week after ADR injection. Bounding box color indicates the condition of the tubule, yellow color indicates normal tubule, green color (labeled with “A”) indicates atrophy tubule, and blue color (labeled with “H”) indicates hypertrophy tubule.

Figure 11 shows the statistical results from the analysis of the complete dataset. The statistics indicates a general trend of increase for average tubule diameter, diameter range, and percentage of hypertrophic tubules. The statistics also shows a fluctuation in these parameters during the midstages (week 2 to 6). We hypothesize it could be due to either variation of different rat’s reaction to ADR injection, sampling bias, or the recovery of the organ before late-stage CKD (week 7 to 8). Previous studies have observed renal function recovery in drug-induced acute kidney injury.⁵³ Tubular density shows an initial increase (week 2 to 4) compared to normal, followed by a decrease in week 5 to 7. This could be explained by the closure of some atrophic tubules.

Fig. 11

The trends of mean tubule diameter, percentage of hypertrophic tubules, tubule diameter range, and tubule density as CKD progression. Matrix of $t$-test between each week is attached to the upper right corner of each plot. The dark regions indicate the $t$-tests between two weeks produce $p$ value above 0.05 (insignificant), and the bright regions indicate $t$-tests between two weeks produce $p$ value below 0.05 (significant).