2.1.

Human Sample Preparation

Freshly resected colon and rectum samples obtained from patients undergoing surgery at Washington University School of Medicine were imaged immediately after surgery. Patients with known benign neoplasia (polyps) as well as malignancies (adenocarcinoma) were eligible for imaging. Cancer patients who had received preoperative treatment with chemotherapy and /or radiation were also included. The study was approved by the Institutional Review Board at Washington University (#201707066). Informed consent was obtained from all patients. Specimens were obtained from the operating room as previously described.⁶

A total of 23 tissue samples were imaged in the pilot study using the PAT/US system. This included untreated colorectal adenocarcinomas ($n=12$), precancerous polyps ($n=6$), colorectal cancer following chemotherapy or radiation and chemotherapy ($n=4$), and postpolypectomy ($n=1$). Two treated patients have achieved complete pathological response and two partial response. The majority of patients underwent hemicolectomy for cancer and were found to have malignancy on histologic analysis (Table 1).

Table 1

Summary of specimens.

2.2.

Co-Registered Ultrasound-Guided Photoacoustic Tomography System

Details of the real-time, co-registered PAT/US system used in this study were discussed previously.^11,22 The system consists of three main parts: a Ti:sapphire laser (Symphotics TII, LS-2134, Symphotics, Camarillo, California) optically pumped with a Q-switched Nd: YAG laser (Symphotics TII, LS-2122), an optimized optical fiber-based light delivery system,²³ and a commercial US system (EC-12R, Alpinion Medical System, Republic of Korea) used for acquiring the corresponding US and PAT data. With this system, pulsed laser light (pulse duration: 10 ns, pulse repetition rate: 15 Hz, $20\mathrm{mJ}/\text{pulse}$ at 750 nm wavelength) was delivered to tissue placed on a two-dimensional motorized scanning stage. B-scan images were acquired for four wavelengths (730, 780, 800, and 830 nm) at each area of imaged tissue among the 23 specimens. The overall scanning area varied from 1 to 3 cm while image acquisition time for a region of interest (ROI) was 12 to 15 s with a frame rate of 15 frames per second. For most of the samples, care was taken to acquire images from abnormal tissues as well as corresponding normal areas, as marked by a trained pathologist.

2.3.

Extraction of Functional, Spectral, and Textural Features

Several functional, spectral, and textural features were extracted from the PAT and US data and images as given in Table 2.

Table 2

Abbreviations.

2.3.2.

Spectral features

Ultrasound images were employed to select a proper ROI corresponding to the lesion for PAT spectral feature calculations.^9,24 First, PAT beam lines with a maximum value close to the background noise level of our co-registered US/PAT system (60 mV) were ignored. The rest of the beam lines were gated by a hamming window, and then their FFT in $-10\text{-}\mathrm{dB}$ frequency range were calculated. Moreover, to cancel the frequency response of the transducer and electrical receiving system,²⁴ the spectra of PAT beams were normalized to the spectra of an approximate point-like target (a $250\text{-}\mu \mathrm{m}$ black thread orthogonal to PAT imaging place with a varied distance to transducer from 0.5 to 7 cm and a step of 0.25 cm).

After calibrating our data, each of the calibrated PAT spectra was fitted linearly. The mean spectral slope (SS), midband fit (MBF), and 0.5-MHz spectral intercept [0.5-MHz SI (PAT)] were then calculated (Fig. 1). We chose 0.5-MHz spectral intercept as a feature instead of 0 MHz because the lower bound of our transducer in PAT mode is $\sim 0.5\mathrm{MHz}$.

Fig. 1

(a), (b) Co-registered rHbT and US images of (a) a cancerous and (b) a normal colon sample. (c), (d) Calibrated PAT power spectra along with their fitted lines in the regions marked with the angular dashed lines.

US spectral features were also calculated. To do this, similar method as PAT spectral features extraction was followed. The only differences were: first, the analysis was performed in the frequency range of 3.5 to 7 MHz, which is the $-10\text{-}\mathrm{dB}$ frequency range of the transducer in US mode. Second, the calibration was performed using a reference gelatin-based phantom constructed in our lab.²⁵

2.3.3.

PAT image features

After visual inspection of PAT frames of malignant and normal colon samples, we noticed that the textures of these images looked different between the two types of samples. To confirm our observation, we extracted PAT image features from available image frames. To do so, a proper ROI was first chosen. To find the center of this ROI, the region surrounding the lesion was determined based on the US image, and the Radon transforms at the two angles of 0 deg and 90 deg of the PAT image in this region were calculated. Each of these Radon transforms was then normalized to its own maximum values and a Gaussian curve was fitted to each of them. The center of the square ROI where the image analysis was performed was determined by the means of these two Gaussian curves, and its size was assumed to be 1 cm for all cases (Fig. 2).

Fig. 2

ROI selection for image analysis. The region covering the lesion is determined (left), and the normalized Radon transforms at 0 and 90 deg are calculated. A Gaussian curve is then fitted to each Radon transform. The center of the square ROI where the image analysis is performed is determined by the means of these two Gaussian curves and its size is 1 cm.

Textural features of the normalized PAT images were calculated in the specified ROI.²⁴ The first step in calculating these features is to construct a gray-level co-occurrence matrix (GLCM).²⁶ GLCM quantifies how the pixels are connected in the image. The size of this matrix was chosen as $16\times 16\text{pixels}$. The value of pixel $(i,j)$ of this matrix was chosen to be the number of times that gray levels $i$ and $j$ are adjacent to each other in the PAT image. Note that we assumed that the two gray levels g1 and g2 are adjacent if g1 is positioned at the immediate left of g2. After constructing the GLCM matrix, four textural features were calculated for each PAT image frame using the following equations:

Eq. (3)

$$\text{contrast}=\sum _{|i-j|=0}^{N-1}{|i-j|}^2\sum _{i=1}^N\sum _{j=1}^{N}c(i,j),$$

Eq. (4)

$$\text{correlation}=\frac{\sum _{i=1}^N\sum _{j=1}^N(i-{\mu }_i)(j-{\mu }_j)c(i,j)}{{\sigma }_i{\sigma }_j},$$

Eq. (5)

$$\text{energy}=\sum _{i=1}^N\sum _{j=1}^{N}c{(i,j)}^2,$$

Eq. (6)

$$\text{homogeneity}=\sum _{i=1}^N\sum _{j=1}^N\frac{c{(i,j)}^2}{1+|i-j|},$$

where $c(i,j)$ is the value of the $(i,j)$ pixel of the GLCM, $N$ is the dimension of this matrix, and $\sigma$ and $\mu$ are the standard deviation and mean for row $i$ or column $j$ of the GLCM.

The standard deviation of the mean radon transform (Sig_rad) was the last image feature that was calculated in this study. To calculate this feature, the radon transform of the non-normalized PAT image at angles 0 deg to 90 deg with a step of 1 deg were calculated and an average was taken over these transforms. Then a Gaussian curve was fitted to this mean radon transform and the standard deviation of this curve was measured.

2.4.

Feature Selection and Classification

A two-step approach was used to select features most likely to differentiate normal from untreated malignant tissue. In the first step, all previously discussed PAT/US features were tested in univariate analysis between untreated cancer and normal regions, and $p$ values were generated by two-sample two-sided Student’s $t$ tests. The features where $p>0.05$—which we concluded a priori not to be significantly associated with malignancy—were excluded from the classification model (Table 3).

Table 3

Significance testing of individual covariates as related to tissue diagnosis.

Next, logistic a general logistic mode (GLM) and a support vector machine (SVM) were used to evaluate the strength of association of each feature with the ultimate tissue diagnosis, and a prediction model was then constructed with significant covariates. In total, 18 areas selected from 18 specimens and 12 malignant areas from 12 untreated cancer specimens were used to construct and evaluate the prediction models. Out of these, 12 normal and 8 malignant areas were used for prediction model derivation and the rest (6 normal and 4 malignant areas) for internal model validation. The receiver operating curve (ROC) and the area under the curve (AUC) were used to evaluate the accuracy of the model. Finally, a second prediction model was constructed without rHbT to determine how limiting the PAT/US device to a single wavelength would affect identification of malignancies.

3.1.

Qualitative Analysis: Baseline Characteristics of US and PAT Images

The colorectal tissues are composed primarily of fluid, lipid, collagen, and muscle. The general architecture (from superficial to deep) in a normal specimen is mucosa (fluid-filled cells surrounded by lipid bilayers), submucosa (largely composed of extracellular collagen matrix and some muscle fibers), muscularis propria (muscle), and adipose tissue (lipid). In malignancy, the individual cell types are similar but the architecture is distorted as cancerous cells of mucosal origin penetrate into the deeper layers of the organ. As these cells invade, the organized structure of the tissue is lost.

Figure 3 shows specimen photographs, US images, co-registered PAT/US rHbT maps as well as histologic images from two representative regions of normal colon samples [Figs. 3(a)–3(d) and 3(e)–3(h)] and two colorectal malignancies [Figs. 3(i)–3(l) and 3(m)–3(p)]. The white arrows indicate the scanning direction along which B-scans were recorded at four different wavelengths (the imaging plane is perpendicular to the scanning direction). In the standalone ultrasound images, the normal layered structure of the colorectal wall is clearly delineated [Figs. 3(b) and 3(f)]. In the presence of malignancy, however, this organized structure is distorted by the tumor and loses the clear delineation of mucosal, submucosal, and muscular layers. These findings mirror the differences in histology among the specimens; in contrast to the ordered layering of the normal colonic wall [images 3(d) and 3(h)], the tumors appear disorganized with destruction of the underlying colonic architecture [images 3(l) and 3(p)].

Fig. 3

Color photograph, US image, rHbT map, and H&E image from representative areas of (a)–(h) two normal regions and (i)–(p) two malignant regions of pretreatment colorectal cancer tissue. Red arrows identify blood vessels within the histologic images.

Additionally, the rHbT maps computed from coregistered PAT/US images of benign regions show significantly lower rHbT signal [Figs. 3(c) and 3(g)] compared to the malignant lesions [3(k) and 3(o)]. As demonstrated in these representative images, normal tissue was found to have almost no detectable rHbT signal. In contrast, malignant tissue showed much higher concentrations of hemoglobin around the tumor bed. Again, these findings appear corroborated by histologic examination. In comparison to the relative paucity of large blood vessels in normal tissue, the malignancies were more vascular and contained large blood vessels [red arrows in images 3(l) and 3(p)].

It is interesting to note that fatty tissues have limited PAT signals in the outer portions of the specimens. This is not surprising since we are specifically targeting hemoglobin—which is not concentrated in fatty tissue—as our chromophore of interest and therefore image within the 730- to 830-nm wavelength range. Additionally, all PAT images are displayed with the same dynamic range of $-10\mathrm{dB}$, so anything below this level is not displayed. The fatty tissue, due to its lack of vascular structures, falls below this range.

3.2.

Evaluation of Treated Tumors

Figure 4 shows corresponding images from a representative untreated colon cancer sample [Figs. 4(a)–4(d)], a colon tumor treated with preoperative chemotherapy [Figs. 4(e)–4(h)], and a rectal malignancy that received radiation and chemotherapy prior to surgical resection [Figs. 4(i)–4(l)]. The untreated tumor displayed findings consistent with other untreated colorectal cancers: loss of layered wall structure, increase in rHbT signal, and increased vascularity throughout the tumor bed. However, treatment appears to reverse these changes. For example, PAT signal appears to diminish with chemotherapy [image 4(g)] and even disappear altogether with complete destruction of the tumor [image 4(k)]. Additionally, ultrasound imaging demonstrated a return to the normal wall structure with complete tumor destruction [image 4(j)]. Histologic comparison among specimens also correlated with these findings; reduction in vasculature along with return to a semiorganized mucosal structure was noted throughout the treated specimens.

Fig. 4

Color photograph, US image, rHbT map, and H&E image from representative areas of (a)–(d) a pretreatment colorectal cancer, (e)–(h) a post-treatment colorectal cancer tissue with residual disease, and (i)–(l) a post-treatment colorectal cancer tissue without residual disease.

3.3.

Quantitative Analysis

In addition to the above qualitative comparisons, a total of 23 areas obtained from 12 untreated malignant tumors, 6 polyps (one has a small invasive component), 2 post-treatment complete responders, 1 no residual tumor cell following prior polypectomy, and 2 post-treatment nonresponders, as well as 18 normal areas from specimens of normal regions were used for quantitative feature extraction. Thus a total of 41 areas were used in Fig. 5. Note that one tumor area was selected from each specimen. Five specimens did not have normal regions or normal regions were too close to the tumor for imaging, therefore, 18 normal areas were selected from 18 specimens. All tumor and normal regions were identified by the attending pathologist.

Fig. 5

Boxplots of (a) total hemoglobin, (b) the mean SS from PAT spectra, (c) 0.5-MHz spectral intercept from PAT spectra, (d) 0.5-MHz spectral intercept from US spectra, (e) energy from the second-order statistics of PAT images, (f) homogeneity from the second-order statistics of PAT images, and (g) standard deviation of the mean radon transform.

Figures 5(a)–5(g) show the boxplots of the seven features calculated from the functional, spectral, and image differences between the different types of colorectal tissue. The $n$ number given in the plots corresponds to the total number of areas. The malignant regions demonstrated elevated rHbT, 0.5-MHz SI (PAT), and 0.5-MHz SI (US) score smaller (less negative) compared to normal and precancerous regions. For SS (PAT), the malignant regions score below normal and precancerous polyps. Treated tumors with complete response were found to have similar scores to normal tissue, whereas treated regions with residual cancer have scores similar to untreated cancers. Due to the limited number of treated cancers, statistics were not performed for these two treated categories.

To distinguish untreated malignant from normal colon tissues, GLM and SVM classifiers were established. These classifiers were developed using the independent features with a $p$-value $<0.05$ between malignant and normal colon tissues. To determine if two features are independent, a spearman’s correlation was calculated between each pair of features (Table 4). To train each classifier, we first used the feature with the lowest $p$-value and then added other features to the feature set one by one. We continued inclusion of the features to the feature set until no increase in the AUC value for the testing data set was observed. We found that when rHbT is included in the feature set, the best performance of both GLM and SVM classifiers (the highest AUC value for the testing data set) is achieved when rHbT and 0.5-MHz SI (PAT) are employed to train the classifier although SS (PAT) has a lower $p$-value than 0.5-MHz SI (PAT). Adding other features did not improve the AUC for the testing data set.

Table 4

The correlation between significant features used in this study.

Figure 6 shows the ROC curves and AUC values of the training (left) and testing (right) data sets using GLM (top) and SVM (bottom) classifiers. As shown in this figure, when the features set include just rHbT, the AUC value for the training and testing data sets are 0.95 and 0.93 for both classifiers, respectively. Adding 0.5-MHz SI (PAT) to the features set, results in a significant improvement in the AUC values for both training and testing data sets (0.97 and 0.95 for the training and testing data sets for both classifiers, respectively). The three image features (Sig_rad, homogeneity, and energy) did not improve the AUC values for both training and testing data sets.

Fig. 6

The ROC curve and their associated AUC values for the training and testing data sets in the presence of rHbT in the feature set. (a), (b) GLM classifier performance. (c), (d) SVM classifier performance.

Finally, the performance of GLM (top) and SVM (bottom) classifiers without rHbT (the single-wavelength model) are presented in Fig. 7. Note that although the difference between some of the PAT image features in malignant and normal samples is statistically significant, none of these features improve the AUC for the testing data sets. The best performance of GLM classifier is achieved when the only spectral feature of SS(PAT) is included in the feature set. The best performance of SVM classifier is achieved when spectral features of SS(PAT), 0.5-MHz SI(PAT), and 0.5-MHz SI(US) are included in the feature set. The testing AUC in this case is 0.89 for the GLM classifier and 0.91 for the SVM classifier.

Fig. 7

The ROC curve and their associated AUC values for the training and testing data sets in the absence of rHbT in the feature set. (a), (b) GLM classifier performance. (c), (d) SVM classifier performance.