*in-vitro*animal experiment to track mucosal tissue inflammation over time, using the area A

_{z}under receiver operating characteristics (ROC) curve as a performance measure. Combination of all the features results in an A

_{z}value up to 1 for the simulated data, and A

_{z}>0.927 for the phantom data. For the tissue data, the best performances for differentiation between pairs of various levels of inflammation are 0.859, 0.983, and 0.999.

## 1.

## Introduction

About 90% of all cancers arise from epithelial cells that cover surfaces of the body and the lining of the internal organs.^{1} Malignancies occurring in the epithelium lining the internal surfaces of organs (e.g., aerodigestive tract or colon) typically are discovered late in the course of the disease, usually because they are “hidden” from the physician, and are not usually found in a routine physical exam. Dysplasia is a precursor tissue change prior to epithelial cancer that can be reversible if it is detected at very early stages. Dysplasia is invisible to the eye, and can only be detected using biopsy.^{1} Investigators are working on alternatives for biopsy to detect precancerous conditions.^{2, 3} Noninvasive and minimally invasive optical techniques are becoming staples of modern medical technology.^{4} Dysplasia can occur less than
$1\phantom{\rule{0.3em}{0ex}}\mathrm{mm}$
in size. For such small lesions, optical techniques are a good choice for detecting the structural changes in the tissue that cannot be easily detected visually.

A variety of optical techniques has been developed and include reflectance,^{5, 6, 7, 8, 9} fluorescence,^{10, 11, 12} light scattering spectroscopy (LSS),^{13, 14, 15, 16} time-resolved autofluorescence,^{12, 17} polarization imaging,^{18, 19, 20} optical coherence tomography (OCT),^{21, 22, 23} diffused optical tomography (DOT),^{24} near-infrared imaging,^{25} Raman scattering spectroscopy,^{26, 27, 28} and frequency domain spectroscopy.^{29} Of these, reflectance and fluorescence methods are of particular interest for their application to early detection of cancer, particularly those of the aerodigestive tract. For example, in the work reported in Ref. 8, cell size and density are determined through backscattered light using spectroscopy. The scattered light is modeled as Mie scattering by surface cell nuclei. In contrast, in Ref. 30, precancerous detection is performed through fluorescence signatures. Since each of these optical modalities has its own points of advantage, combinations of such methods often produce specificity and sensitivities that are unattainable by individual techniques.^{7}

Optical techniques do not necessitate tissue removal and analysis can be made in real time, hence offering the potential of detecting mucosal changes at the microstructural, biochemical, and molecular level. The detailed spatial intensity distribution of light scattered by an individual particle is a complex function of the particle’s size, shape, and orientation with respect to the wavelength as well as the incident illumination.^{31} Among other important measures, scattered light provides an objective measure of epithelial nuclear enlargement and crowding, which are the most significant characteristics of dysplasia and early cancer.^{1}

Our goal is to establish a textural model that models the reflected scattered signal and shows its sensitivity to differential changes in the morphology and physical characteristics of the tissue being optically scanned. Textural models have been successfully used by us to model radio frequency (rf) ultrasound images of breast tissue,^{32, 33, 34} where we have proven the model to be quite effective for discrimination between benign, normal, and malignant breast tissues (the
${A}_{z}$
values were in the range of 0.862 to 0.999). Similar to our work on rf modeling of breast tissue, we adopt a stochastic decomposition method (SDM) to model the scattered signal optically reflected from mucosal tissues to track down differential changes in the morphology and physical characteristics of the imaged (scanned) tissue. One of the original contributions of this work is the application of the SDM model to optical data, which is totally new to this application field. A second contribution is the systematic analysis and correlation between the signal parameters in predicting and tracking morphological and physical characteristics of epithelial tissue as it undergoes changes. Simulations and phantom data are used to test the validity of the model. The test scenarios are driven with biological applications (such as dysplasia) in mind. Differentiation based on nuclei sizes and densities (a trademark of nuclear crowding in dysplasia) were embedded in the choice of the simulation parameters of the light scattered from the cells mimicking dysplasia based on the Monte Carlo (stochastic) simulation of Mie scattering. Some of the parameters in our stochastic Mie scattering simulator relate directly to the tissue morphology, such as: number of cells, size of cells, and location of cells; while others relate to the physical staining characteristics of the tissue, which affect and are affected by the absorption and refraction indices. The model is also tested using phantoms with nuclei sizes and density scatterers (spheres) varying in accordance with biological epithelial tissue as it undergoes dysplastic changes. We used phantom data to see the performance of the technique to discriminate between two different morphologies reflecting tissue size of cell nuclei sizes mimicking normal and dysplastic epithelium. The simulations and phantom experiments show that some of the SDM parameters are extremely effective in picking up changes in tissue morphologies.

To further illustrate the feasibility of the model, we have considered a preliminary animal model from mouse colon *in vitro* with different levels of induced inflammation, where we test the ability of the SDM to differentiate between the various stages of inflammation as it develops in the mouse colon epithelium. This is of particular relevance, as there is a functional relationship between inflammation and cancer.^{35} Balkwill and Mantovani^{36} described that issue as: if genetic damage is the “match that lights the fire” of cancer, some types of inflammation may provide the “fuel that feeds the flames.”^{36} For instance, ulcerative colitis (UC) associated colorectal cancer follows a histological sequence starting in the inflamed mucosa as a hyperplastic lesion develops through dysplasia into adenocarcinoma.^{37, 38} This is sometimes summarized as the “inflammation-dysplasia-carcinoma” sequence.^{39}

## 2.

## Stochastic Decomposition Method and the Extracted Features

The decomposition of the LSS image (2-D scan) is consistent with the general decomposition of regular stochastic fields into predictable (the specular field) components and unpredictable (the diffused field) components, known in the literature of stochastic processes as the World decomposition, used to decompose the scattered or reflected signal into its two components.^{33, 40} The specular scattering component is modeled by periodic or quasiperiodic image field of point scatterers corresponding to the cell boundaries, whereas the diffuse component is modeled by an autoregressive field, which corresponds to a linear filter driven by white noise:

The resolvable scattering structure $c\left(n\right)$ can be viewed as a summation of smeared modulated delta functions located at the resolvable scatterers’ locations, and of random strength, the coherent component $c\left(n\right)$ can be modeled by a superposition of Gaussian modulated sinusoids:

## Eq. 3

$$c\left(n\right)=\sum _{l=0}^{{N}_{e}}\frac{{A}_{l}}{{\sigma}_{cl}\sqrt{2\pi}}\phantom{\rule{0.2em}{0ex}}\mathrm{exp}\left\{\frac{{(n-{m}_{i})}^{2}}{2{\sigma}_{cl}^{2}}\right\}\mathrm{cos}\left({w}_{c}n\right),$$The SDM first checks for the existence of coherent scatterers by testing the hypothesis of Rayleigh scattering, which can be achieved by using the nonparametric Kolmogorov-Smirnov (KS) test for color field.^{41} The rejection of the Rayleigh scattering hypothesis indicates the existence of a coherent component. If a coherent component exists, the signal is decomposed into its coherent
$c\left(n\right)$
and diffuse
$d\left(n\right)$
components using the wavelet decomposition described in Ref. 32. The SDM consists of taking the wavelet transform of the scattered signal and thresholding its energy. The decomposition uses the continuous wavelet transform and was thoroughly described and tested on simulated ultrasound rf data.^{32, 33} After the decomposition, we fit the autoregressive (AR) model to the diffuse part and extract the parameters of the diffuse part as well as the coherent part. If the test of presence of a coherent component fails, the signal is totally diffused and is modeled as an AR process.

The parameters for the coherent component are as follows.

*Number of detected coherent scatterers N _{c} in a region of interest (ROI).*The locations of large fluctuation of the scale-averaged wavelet power around its mean value are the locations of the coherent scatterers in the probed region under study.

^{33, 34}The number of the fluctuations that exceed a threshold is the number of coherent scatterers in the tissue. The number of coherent scatterers is related to the number of resolvable scatterers in the tissue (single scatterers). This parameter is correlated with the size of the window. The larger the size of the window is, the more coherent scatterers might be detected.

*Mean energy of the coherent scatterers
$E$
*. This parameter reflects differences in energy due to changes in index of refraction that affect the absorption and reflected energy.

The parameters and features of the diffuse component are estimated from the model parameters and are given as follows.

*Residual error variance of the diffuse component*
${\sigma}^{2}$
. The equivalent of the energy of the coherent scatterers for the diffuse component is the unconditional variance. The unconditional variance is related to the residual error variance (conditional variance of the diffuse component) through the coefficients of the autoregressive model (AR parameters). The values of the unconditional variance are reflected mainly on the residual error variance.^{33, 34} They are also weakly reflected on the AR parameters. It follows that the same argument made for the mean energy of the coherent scatterers applies for the residual error variance as well, because if the resolution increases drastically, the previously unresolved scatterers will be resolvable. In our previous research for breast cancer characterization using ultrasound, this parameter showed excellent discrimination ability between benign, normal, and malignant tissues.^{33}

*Rayleigh scattering degree of the diffuse component*
$D$
. This parameter describes the discrepancy between the empirical distribution of the innovation process of the diffuse component and the Gaussian distribution by calculating the KS distance. The KS distance is the distance between two cumulative probability distribution functions; it is the distance used to test the hypothesis that given data follow a certain nominal distribution. The lower the value of
$D$
, the closer to Rayleigh scattering the diffuse component is.^{33, 34} In our previous research with ultrasound backscattering, this parameter was found to be tightly connected with the density of scatterers per resolution cell; large scatterers per resolution cell caused Rayleigh scattering and less scatterers caused deviation from that.^{33} Based on our past experience, we anticipate this parameter to be able to differentiate between the various stages of dyplasia formation as it directly relates to formation of new dysplastic cells and the increase/decrease in the dysplastic cell density that directly reflects on that KS distance.

*Normalized correlation coefficient of the diffuse component*
${\rho}_{N}$
This parameter reflects the correlations between neighboring diffuse scatterers in the tissue, and depends on the density of the diffuse scatterers in the tissue in relation to the sampling period. The normalized correlation coefficient is the determinant of the covariance matrix normalized to set the diagonal of the matrix equal to 1.^{33} The determinant of the normalized covariance matrix is real, nonnegative, and equals 0 if the data are linearly dependent, and equals 1 if the data are independent.

The parameters can be analyzed either individually and/or in combination to construct a detailed map of a tissue region. When using parameters in combination, global decision fusion must be used.^{34} The global decision fusion we use is explained in detail in the following paragraphs.

Let
${Y}_{i}$
be the
$i$
’th feature,
${H}_{0}$
be the null hypothesis (normal tissue), and
${H}_{1}$
be the alternative hypothesis (cancerous tissue). The local Neyman Pearson (NP) decision rule with the highest detection rate among all other decision rules for any fixed false alarm rate is^{34}

## Eq. 4

$$\frac{p({Y}_{i}\mid {H}_{1})}{p({Y}_{i}\mid {H}_{0})}\underset{{H}_{0}}{\overset{{H}_{1}}{\gtrless}}{\zeta}_{i}.$$## Eq. 5

$${X}_{i}=1\phantom{\rule{0.3em}{0ex}}\mathrm{iff}\phantom{\rule{0.3em}{0ex}}{T}_{i}>{\zeta}_{i}\phantom{\rule{1em}{0ex}}\text{and}\phantom{\rule{1em}{0ex}}{X}_{i}=0\phantom{\rule{0.3em}{0ex}}\mathrm{iff}\phantom{\rule{0.3em}{0ex}}{T}_{i}<{\zeta}_{i},$$^{34}

## Eq. 6

$${T}_{i}=\frac{{({Y}_{i}-{\mu}_{0})}^{2}}{{\sigma}_{0}^{2}}-\frac{{({Y}_{i}-{\mu}_{1})}^{2}}{{\sigma}_{1}^{2}},$$^{34}

## Eq. 7

$$\frac{p({\stackrel{\u20d7}{Y}}_{H}\mid {H}_{1})}{p({\stackrel{\u20d7}{Y}}_{H}\mid {H}_{o})}\underset{{H}_{0}}{\overset{{H}_{1}}{\gtrless}}\zeta ,$$^{34}:

## Eq. 8

$${T}_{g}={(\stackrel{\u20d7}{Y}-{\stackrel{\u20d7}{\mu}}_{0})}^{T}{\Sigma}_{0}^{-1}(\stackrel{\u20d7}{Y}-{\stackrel{\u20d7}{\mu}}_{0})-{(\stackrel{\u20d7}{Y}-{\stackrel{\u20d7}{\mu}}_{1})}^{T}{\Sigma}_{1}^{-1}(\stackrel{\u20d7}{Y}-{\stackrel{\u20d7}{\mu}}_{1}).$$Since the NP has the highest probability of detection for a fixed false alarm rate among all other decision rules, it therefore follows that the power (detection rate) of the global likelihood ratio test
$(1-\beta )\u2a7e(1-{\beta}_{i})$
,
$i\u220aH$
. Hence, the global receiver operating characteristics (ROC) curve, which is a plot of probability of detection versus probability of false alarm, bounds from above all the local ROC curves. Simply stated, this proves that the addition of further information (sum of features) improves the probability of detection over using the individual features alone, and that is true for any given false alarm rate.^{34} This means that the fusion of the various SDM features will result in better detection than either alone.

We use the area under the ROC curve
${A}_{z}$
as a metric, since it integrates out the dependence of the performance on which false alarm rate to use. A sample ROC curve and sample class distributions from which the ROC curve is determined are shown in Fig. 1
.
${A}_{z}$
values close to 0.5 indicate very weak discrimination power between two classes (total chance decision rule amounting to coin tossing), whereas values close to 1 indicate very strong discrimination power and yield very high sensitivity and specificity values.^{34}

## 3.

## Details of Experiments and Data Collection

In the experiments either based on simulations or on phantom data, we select the parameters and their values to mimic the interrogation of a biological epithelium tissue that is undergoing neoplastic changes, either dysplastic or of a malignant nature. Dysplasia is a term that refers to a precancerous condition. Malignant neoplastic changes mostly follow preexisting dysplastic changes. Removal of adverse environmental stimulus leads to restoration of normal cell growth patterns, which makes dysplasia reversible if detected in its very early stages. Dysplasia is recognized by cells altering their appearance (cytology). As tissue becomes dysplastic, the nuclei enlarge and become crowded.^{6} Healthy tissue epithelial nuclei have a characteristic diameter of
$4\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}7\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
and are arranged in neat rows.^{8} In dysplastic epithelium, the cells proliferate and the nuclei enlarge and appear darker when stained.^{1} Nuclei can be as large as
$20\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
in diameter.^{8} Normal intestinal cells are characterized by uniform nuclear size distribution where malignant cells have larger nuclei and more variation in nuclear size.^{18} In Ref. 18, for a normal tissue sample, the average diameter was found to be
$4.8\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
, a standard deviation of the sizes was
$0.4\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
. For the cancerous tissue sample, the corresponding values were 9.75 and
$1.5\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
. The detection of dysplasia is a challenging problem, since it does not form macroscopically observable structures that can be picked during colonoscopy.

Finally, to further show the feasibility of the model, we have also considered a preliminary *in vitro* animal model of a mouse colon with different levels of induced inflammation. Here we test the ability of the SDM to differentiate between the various stages of inflammation as it develops in the mouse colon epithelium.

## 3.1.

### Simulation of Scattered Light from Cluster of Spheres Subjected to Changes in Morphology and Physical Characteristics

The simulations are based on linear arrays of light scattering data at a given wavelength. The results on simulated overlapping and nonoverlapping 1-D scans (256 points) of data obeying Mie scattering have been obtained. The theoretical details for the simulations are given in the Appendix (Sec. 6). Each intensity point [See Eq. 14] on the 1-D scan results from a cluster of $N$ multiscattering spheres located at given positions $(x,y,z)$ within a cube extending from $-20\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}+20\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ on all three axes. Unless otherwise stated, the term “resolution cell” refers to this volume. Since we do not deal with overlapping spheres, the sphere locations are selected to be sufficiently spaced. The spheres were of a given size and index of refraction (RI) Figure 2 shows an ensemble of sphere scattering clusters, with each cluster comprising $N$ spheres, each with diameter $d$ . The strength of the incident light $I$ at a fixed wavelength $\lambda $ ( $580\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ in the experiment) and the direction of the detector $(\theta ,\phi )$ were control parameters in the simulator. The scattering matrix corresponding to a direction of (45 and $90\phantom{\rule{0.3em}{0ex}}\mathrm{deg}$ ) is calculated. The output of the simulator is the scattered light received at the detector from the ensemble of the scattering spheres obtained using the Mueller scattering matrix [see Eq. 11]. The scattered light was calculated at various positions of the resolution cell (centroids of the cube). The contour map is rectangular, with $\theta $ and $\phi $ values representing the vertical and horizontal coordinates. The output of the program is the Mueller scattering matrix map written sequentially for each set of $\theta $ and $\phi $ values. The simulator calculates the scattering matrix for $nt+1$ values of $\theta $ between 0 and $\pi $ , and $nt+1$ values of $\phi $ between 0 and $2\pi $ . The output considered out of this $nt+1$ by $nt+1$ matrix always corresponds to $45\phantom{\rule{0.3em}{0ex}}\mathrm{deg}$ ( $\theta $ value) and $90\phantom{\rule{0.3em}{0ex}}\mathrm{deg}$ ( $\phi $ value). The output of the simulation gives the elements of the Mueller scattering matrix of the scattered light, corresponding to $\theta =45\phantom{\rule{0.3em}{0ex}}\mathrm{deg}$ and $\phi =90\phantom{\rule{0.3em}{0ex}}\mathrm{deg}$ [Eqs. 12, 13]. The intensity is calculated from the scattering matrix using Eq. 14.

To introduce randomness into Mie scattering,^{16} we slightly and randomly perturb the spheres’ position coordinates, and number and refractive index (RI) of the scattering spheres around nominal values defining a given tissue structure. The coordinates of the first cluster ensemble are specified and the coordinates for the remaining 256 clusters are calculated, generating from the first cluster by adding a random position component to each of the
$N$
spheres in the first cluster using a Gaussian random number generator. The position of each sphere in the first cluster using a Gaussian random number generator. The position of each sphere in the new cluster is computed as follows:

## Eq. 9

$${(x,y,z)}_{i}={(x,y,z)}_{1}+\text{random}\phantom{\rule{0.3em}{0ex}}\text{coordinates}\phantom{\rule{0.3em}{0ex}}\text{from}\phantom{\rule{0.3em}{0ex}}\text{random}\phantom{\rule{0.3em}{0ex}}\text{number}\phantom{\rule{0.3em}{0ex}}\text{generator}.$$We have considered resolution cells and we picked parameters [cell sizes, number, and index of refraction (RI)] that are close in values to normal and crowding epithelial nuclei that are encountered in dysplasia. Two types of cells are simulated, normal cells and diseased cells. For the simulation of the *normal cells*, the cluster volume is assumed to contain 9 to 11 cells of
$5\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
diameter. For the simulation of the *diseased cells*, each cluster ensemble is considered to consist of four cells of
$12\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
diameter. The index of refraction (RI) of the tissue is generally in the range of 1.33 to 1.5.^{42, 43} RI takes a minimal value when the content of water in tissue is maximal and vice versa.^{42, 44} We used two different indices of refraction,
$R{I}_{1}$
and
$R{I}_{2}$
(
$R{I}_{1}=1.38+0.005i$
and
$R{I}_{2}=1.41+0.010i$
), values that are consistent with many other studies in the literature.
^{42, 45, 46, 47, 48, 49, 50, 51, 52, 53}

Three different sets of simulations are considered. In simulation sets 1 and 2, four types of data are created, L1, L2, S1, and S2. L1 stands for large spheres with cell diameter of $12\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ with index of refraction. $1.38+0.005i$ , whereas L2 stands for the same configuration except that the index of refraction is increased to $1.41+0.01i$ . S1 stands for small spheres with cell diameters of $5\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ with index of refraction $1.38+0.005i$ , whereas S2 stands for the same configuration except that the index of refraction is increased to $1.41+0.01i$ . L1 and L2 represent diseased cells, while S1 and S2 represent normal cells. Descriptions for the four types of data are summarized in Table 1 .

## Table 1

Parameter values for simulation sets 1 and 2.

Cell Diameter | Number of Cells | Index of Refraction | ||
---|---|---|---|---|

Diseasedcells | L1 | $12\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ | 4 | $1.38+0.005i$ |

L2 | $12\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ | 4 | $1.41+0.010i$ | |

Normal cells | S1 | $5\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ | 9 to 11 | $1.38+0.005i$ |

S2 | $5\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ | 9 to 11 | $1.41+0.010i$ |

We repeat the experiment (simulation set 2), but now we vary the resolution so that some of the scatterers behave as coherent scatterers as well as diffuse. For simulation set 1, there is no overlap between consecutive samples. However, for this new case, we take samples that are measured every $10\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ , thus producing an overlap between the 256 samples as opposed to the simulation set 1.

Finally, in simulation set 3 we consider the case of one single scatterer per resolution cell as opposed to many (multiscattering), and there is no overlap between consecutive samples. Three types of data are created, A1, A2, and A3. For A1 and A2, there is a single sphere with diameters 30 and $10\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ , respectively, in the resolution cell. All 256 samples are randomized through small random uniform perturbations on the index of refraction around $1.38+0.005i$ : positions of the spheres are fixed from cluster to cluster. A3 has ten spheres (multiscattering case) of $5\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ diameter and index of refraction $1.38+0.005i$ . The descriptions of these three types of data are given in Table 2 .

## Table 2

Parameter values for simulation set 3.

Cell diameter | Number of cells | Index of refraction | ||
---|---|---|---|---|

Single spherescattering | A1 | $30\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ | 1 | Perturbed around $1.38+0.005i$ |

A2 | $10\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ | 1 | Perturbed around $1.38+0.005i$ | |

MultiSpherescattering | A3 | $5\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ | 10 | $1.38+0.005i$ |

## 3.2.

### Details of Phantom Experiments

(Probe, spectrometer, and light source are provided by Drexel Photonics Laboratory members Vitol, Kurzweg and Nabet). Phantom and *in-vitro* animal data are collected using a probe, spectrometer, and light source.^{54} The sample is illuminated with white light, and the reflected light spectrum is collected at equal intervals during scanning using a staging station. Intensities collected at one specific wavelength at each data point from an image (2-D scan) of the sample. Various 2-D scans are formed by the intensities at different specific wavelengths. Our method images the tissue area under investigation, and extracts from the data features (SDM parameters) that change with changes in the tissue morphology. A distinctive feature of our technique compared to DOT, OCT, and LSS is the formation of an image (2-D scan), rather than a spot, from which pertinent data are extracted. In this fashion, obtaining more data cannot only lead to a more confident calculation of the relevant features, such as nuclear size distribution, but can also lead to additional information embedded in the spatial texture that our decomposition technique arrives at by modeling the hidden correlations that are obtained only by interrogating a wide sample area.

In our first set of experiments, we use polystyrene latex microspheres (Polysciences Incorporated, Washington, Pennsylvania) suspended in deionized water.^{54} We used phantom data to see the performance of the technique to discriminate between two different morphologies reflecting tissue cell nuclei sizes mimicking normal and dysplastic epithelium. The microspheres exist in a 2.5% aqueous suspension with a refractive index of 1.59 to 1.60 at
$589\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
. We examined microspheres at diameter sizes of 3 and
$10\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
. A schematic diagram for the experimental system setup is shown in Fig. 3
, and a picture of the actual system (probe, spectrometer, and
$xyz$
stage) is shown in Fig. 4
. We use a bifurcated reflectance probe (Ocean Optics, Dunedin, Florida) consisting of one central fiber used for light delivery and six surrounding fibers used for light collection, each with a core of
$200\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
.^{54} The probe has a special
$30\phantom{\rule{0.3em}{0ex}}\mathrm{deg}$
window at the end to reduce specular reflectance from the sample surface. The probe is connected to the light source (tungsten halogen lamp; LS-1, Ocean Optics, Dunedin, Floria) and the spectrometer (HR4000, Ocean Optics, Dunedin, Florida) having an output of 3648 points with a wavelength step size of
$0.12\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$
. To create the image (2-D scan), the fiber is placed on an
$xyz$
stage, which allows the fiber to be moved and data to be collected at different points on the sample. At each data point, the whole spectrum is recorded. After the spectrum is collected at all the data points, the intensity image (2-D scan) is reconstructed by examining one single wavelength from which key parameters are extracted for each performance evaluation.

We use $3\text{-}\mu \mathrm{m}$ spheres to mimic the normal cells, and $10\text{-}\mu \mathrm{m}$ spheres to mimic the dysplastic cells. The surface area over which the data is collected is $4500\times 4500\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ . Three 1-D scans are taken for each phantom. The scans are taken at 0, 2250, and $3500\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ (the $y$ direction). For each 1-D scan, 51 points are taken $90\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$ apart in the $x$ direction. At each point the whole spectrum (at various wavelengths) is computed. Smoothing is used as a preprocessing step for each point. After the spectrum is smoothed, the intensity values for specific wavelengths are taken. As the phantom was limited in size, order permutations are carried out to generate a set of 32 independent 1-D scans from each of the three 1-D scans, resulting into a total of 99 1-D scans used for classification.

## 3.3.

### Details of Tissue Experiments

The very preliminary animal study is based on the use of mid-colon obtained from three mice fed with 4% dextrin sulfate solution (DSS). The DSS model has been used to study colorectal cancer (CRC). The DSS colitis model can be used in mouse as a model for studying the sequence of chronic inflammation dysplasia in human inflammatory bowel disease.^{55} For normal ulcerative colitis, mice are fed with DSS for one cycle; for chronic inflammation one cycle of DSS is followed by two cycles of normal water. For dysplasia, mice are kept on DSS for one cycle, two cycles of water, again DSS, then followed by water. In our model, what we were measuring was uncontrolled acute inflammation. Prior histological studies have shown that during the first seven days, of DSS feeding, these animals show gradual crypt dropout and acute to mixed inflammation starting on day 3 of DSS feeding.^{56} Colon was excised on days 4 and 7 of DSS feeding. Clinically, after four days, the disease activity index was 1.3, which increased to 2.7 on day 7. In our experiments, the inflammation certainly increased from day to day while on DSS feed. We named the samples as mouse1, mouse2, and mouse3. Mouse1 represents normal mouse colon *in vitro*, mouse2 represents increased inflammation with the presence of increased neutrophils (polymorph nucleus) on day 4, and mouse3 represents further increase in inflammation with the possible presence of mucosal erosions on day 7. The excised colon was mounted in a petri dish containing phosphate buffered saline. The colon was imaged within thirty minutes after mounting. It takes about
$7\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}8\phantom{\rule{0.3em}{0ex}}\mathrm{min}$
to complete the scan. A schematic of the animal study is shown in Fig. 5
.

Three 1-D scans of the entire piece of the colon
$(\mathrm{10,000}\times 3000\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m})$
were taken in a zigzag fashion, using a micrometer stage within
$30\phantom{\rule{0.3em}{0ex}}\mathrm{min}$
after tissue mounting. The data collection process followed a similar protocol as with the phantom data collection. A 3-D figure of the experimental setup for *in-vitro* data collection is shown in Fig. 3. Three scans are taken as in phantom data. The scans are taken at 500, 1500, and
$2500\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
. Each scan consists of 50 points. The step size is taken as
$200\phantom{\rule{0.3em}{0ex}}\mu \mathrm{m}$
. After the data are collected, the same processing steps as with the phantom data are applied.

## 4.

## Results and Discussion

The metric that we used for representing the results is the area ${A}_{z}$ under the ROC curve (plot of probability of detection versus probability of false alarm). As explained in Sec. 2, we use the ${A}_{z}$ value under the ROC curve to determine the performance of our classification system.

## 4.1.

### Performance Evaluation on Simulations

Using the SDM, we calculated five features ${N}_{c}$ , $D$ , ${\sigma}^{2}$ , ${\rho}_{N}$ , and $E$ for each realization. The classification performance of each feature is evaluated separately. For simulation set 1, there is no overlap between consecutive samples. At the resolution used, the test for the existence of coherent scatterers was rejected as expected, since for this case only diffuse scatterers were present. The performances ( ${A}_{z}$ values) of individual diffuse component features for the four cases for simulation set 1 are given in Table 3 . Each value in the table presents the area under the ROC curve for classification between pairs of data. Using a combination of features, the global detector yielded a perfect performance with overall ${A}_{z}$ values greater than 0.991. The combination of performance parameters is in accord with the joint likelihood ratio statistics (NP statistics), which insures the highest detection rate for a fixed alarm rate over any other classifier including the NP based on the individual parameters in isolation (see Sec. 2).

## Table 3

Performance of the system for the extracted features in terms of Az (area under the ROC curve) values for differentiation between the pairs in simulation set 1.

L1-S1 | L2-S2 | L1-L2 | S1-S2 | |
---|---|---|---|---|

$D$ | 0.971 | 0.996 | 0.672 | 0.790 |

${\sigma}^{2}$ | 1.000 | 1.000 | 0.962 | 1.000 |

${\rho}_{N}$ | 0.586 | 0.583 | 0.537 | 0.505 |

All | 1.000 | 1.000 | 0.991 | 1.000 |

It is interesting to note here that the KS distance parameter $D$ , which directly relates to the number of scatterers (spheres), was able to differentiate between structures with different numbers and sizes of scattering spheres but failed (as expected) for the case with same density (L1-L2 and S1-S2) but with different indexes of refraction. On the other hand, since the energy of the scattered light (and hence the residual error variance) depends both on the number of scatterers and the index of refraction, as expected, it performed well in all four cases. Finally, as expected, the normalized correlation coefficient was not able to discriminate between the various cases, as the 256 points resulted from nonoverlapping spheres as we moved along the simulated structure. The 3-D plots of the extracted features are given in Fig. 6 to show the good separation between different pairs of classes in simulation set 1.

For simulation set 2, the performances ( ${A}_{z}$ values) of individual diffuse component features for the four cases are given in Table 4 . Using a combination of features, the global detector yields a perfect performance with overall ${A}_{z}$ values greater than 0.999 for the cases. Unlike the previous case, since there is an overlap between consecutive samples, a coherent component is detected and we have coherent as well as diffuse parameters. In addition, the normalized correlation coefficient becomes a discriminating parameter for cases L1-S1 and L2-S2, where the data between pairs of classes correspond to different numbers of scatterers and size.

## Table 4

Performance of the system for the extracted features in terms of Az values for differentiation between the pairs in simulation set 2.

L1-S1 | L2-S2 | L1-L2 | S1-S2 | |
---|---|---|---|---|

${N}_{c}$ | 0.993 | 1.000 | 0.639 | 0.507 |

$E$ | 0.999 | 1.000 | 0.889 | 0.965 |

$D$ | 0.987 | 0.968 | 0.568 | 0.546 |

${\sigma}^{2}$ | 1.000 | 1.000 | 0.948 | 0.983 |

${\rho}_{N}$ | 0.890 | 0.807 | 0.679 | 0.607 |

All | 1.000 | 1.000 | 0.999 | 1.000 |

To discover the differentiation performance of the method for single and multisphere scattering cases, simulation set 3 is used, and the following results are obtained. The performances ( ${A}_{z}$ values) of the individual features for the pairs A1-A2, A1-A3, and A2-A3 are given in Table 5 . Using a combination of features, the global detector yields a perfect performance with an overall ${A}_{z}$ value of 1 for the cases. For the case of single scattering, as expected, the number of coherent scatterers as well as the KS distance cannot discriminate between cases A1 and A2, since they both have one scatterer per resolution cell. However, they do successfully differentiate when the number of scattering spheres is different (cases A1-A3 and A2-A3). As expected, the parameters (the energy of the coherent and diffuse scatterers) successfully discriminate between all cases, since they depend on both the index of refraction and the number of scatterers per resolution cell. Finally, as anticipated, the normalized correlation coefficient is not able to discriminate between the various cases due to data permutations that destroyed any discrimination based on correlation. The 3-D plots of the extracted features are given in Fig. 7 to show how good the separation was between different pairs of classes in simulation set 3.

## Table 5

Performance of the system for the extracted features in terms of Az values for differentiation between the pairs in simulation set 3.

A1-A3 | A2-A3 | A1-A2 | |
---|---|---|---|

${N}_{c}$ | 0.998 | 1.000 | 0.683 |

$E$ | 0.989 | 1.000 | 1.000 |

$D$ | 0.990 | 0.975 | 0.592 |

${\sigma}^{2}$ | 1.000 | 1.000 | 1.000 |

${\rho}_{N}$ | 0.695 | 0.676 | 0.578 |

All | 1.000 | 1.000 | 1.000 |

## 4.2.

### Performance Evaluation on Phantom

The performance of the system is calculated for each specific wavelength. For each 1-D scan at a given wavelength, the decomposition parameters were computed, hence a set of 99 imaging parameter datasets for the $3\text{-}\mu \mathrm{m}$ spheres and 99 for the $10\text{-}\mu \mathrm{m}$ spheres. The classification is reported to be the best for wavelengths between $500\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}600\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ . The classification performance between the normal versus the dysplastic mimicking cells is evaluated, as in the simulation case, using the area under the ROC curve as the metric for performance. This is shown in Tables 6, 7 for a diffuse scatterer only model versus a diffuse plus coherent scatterer model. The last row in both tables shows, the results of classification by fusing all the imaging parameter sets together. The performance varied from 0.927 to 0.994, which is consistent with the simulation results, and is a very strong result in discriminating between normal-mimicking cells from dysplastic-mimicking ones. As with the simulation results, the best feature is still the residual error variance. The overall performance is above 0.9 for both the simulated and phantom data. It is shown that the performances of the individual parameters remain consistent with the simulation results, and the overall performance of the system is still extremely high $({A}_{z}>0.927)$ .

## Table 6

Performance of the system for the extracted features in terms of Az values for diffuse only model for phantom data at different wavelengths.

500.04nm | 550.05nm | 571.67nm | 650.00nm | |
---|---|---|---|---|

$D$ | 0.694 | 0.551 | 0.580 | 0.519 |

${\sigma}^{2}$ | 0.981 | 0.979 | 0.982 | 0.924 |

${\rho}_{N}$ | 0.500 | 0.576 | 0.561 | 0.504 |

All | 0.990 | 0.982 | 0.989 | 0.927 |

## Table 7

Performance of the system for the extracted features in terms of Az values for coherent+diffuse model for phantom data at different wavelengths.

500.04nm | 550.05nm | 571.67nm | 650.00nm | |
---|---|---|---|---|

${N}_{c}$ | 0.551 | 0.587 | 0.517 | 0.575 |

$E$ | 0.986 | 0.967 | 0.979 | 0.865 |

$D$ | 0.733 | 0.530 | 0.718 | 0.616 |

${\sigma}^{2}$ | 0.984 | 0.979 | 0.982 | 0.924 |

${\rho}_{N}$ | 0.572 | 0.538 | 0.508 | 0.510 |

All | 0.994 | 0.988 | 0.991 | 0.931 |

## 4.3.

### Performance Evaluation on Tissue

The classification is reported for wavelengths between $450\phantom{\rule{0.3em}{0ex}}\text{to}\phantom{\rule{0.3em}{0ex}}650\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ in Tables 8, 9, 10 . As with the phantom results, we observe that the mean energy of the coherent and diffuse scatterer increased as the animals became sicker. The performance is shown in the Tables 8, 9 in terms of areas under the ROC curves ( ${A}_{z}$ values) using again the permutation method used with the phantom data to increase the sample size (i.e., the number of A scans). The best performance using one wavelength $\left(625.24\phantom{\rule{0.3em}{0ex}}\mathrm{nm}\right)$ was at 0.859 combined ${A}_{z}$ value (i.e., using all features) for differentiating between normal tissue and $4\text{-}\text{day}$ inflammation, whereas it was at 0.999 combined ${A}_{z}$ value for differentiating between normal tissue versus day-7 inflammation at the same wavelength of $625.24\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ . We have also run the classification algorithm to see if we can differentiate between the $4\text{-}\text{day}$ inflammations versus the $7\text{-}\text{day}$ ones. The results are shown in Table 10. The best performance using one wavelength was at 0.987 combined ${A}_{z}$ value for differentiating between $4\text{-}\text{day}$ inflammation tissue versus $7\text{-}\text{day}$ inflammation tissue at a wavelength of $575.02\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ . The best performance of the algorithm was obtained at the wavelengths around the peak values of the spectrum. For the best two performing parameters $E$ and ${\sigma}^{2}$ , the wavelengths that show the best detection performances show consistency with each other for all the classification pairs. Also, they show consistency with the corresponding wavelengths for the overall performance (performance generated by fusing all the features). Although we have only considered three mice, we have looked at various area or points within each of these samples (three scans per sample at different locations with 50 points each). In Fig. 8 , we show the scattered data for the best two features for the three mice at $575.02\phantom{\rule{0.3em}{0ex}}\mathrm{nm}$ . The ROC construction and the ${A}_{z}$ values are based on the scattered data points. We can easily see from the scattered data how well separated the data are between the normal and $7\text{-}\text{day}$ inflammation, and how it is less so between the normal and the $4\text{-}\text{day}$ inflammation, leading to high ${A}_{z}$ value under the ROC curve for the former case and lower for the latter. The separation between the $4\text{-}\text{day}$ and the $7\text{-}\text{day}$ inflammations is also quite pronounced, leading to high ${A}_{z}$ value under the ROC curve. The scatterer plot shows the direct proportionality of the features to the disease formation. As the disease level increases, the values of the features’ “mean energy of coherent scatterers” and “residual error variance” increase as well. This result is promising and consistent with the findings of the simulations and phantom data.

## Table 8

Performance of the system for the extracted features in terms of Az values for tissue data (mouse colon in vitro) for classification between case Mouse 1 versus Mouse 2.

500.33nm | 550.02nm | 575.02nm | 625.24nm | |
---|---|---|---|---|

${N}_{c}$ | 0.568 | 0.559 | 0.590 | 0.559 |

$E$ | 0.811 | 0.741 | 0.743 | 0.843 |

$D$ | 0.515 | 0.505 | 0.563 | 0.565 |

${\sigma}^{2}$ | 0.809 | 0.703 | 0.718 | 0.821 |

${\rho}_{N}$ | 0.510 | 0.561 | 0.506 | 0.513 |

All | 0.812 | 0.770 | 0.798 | 0.859 |

## Table 9

Performance of the system for the extracted features in terms of Az values for tissue data (mouse colon in vitro) for classification between case Mouse 1 versus Mouse 3.

500.33nm | 550.02nm | 575.02nm | 625.24nm | |
---|---|---|---|---|

${N}_{c}$ | 0.541 | 0.511 | 0.549 | 0.508 |

$E$ | 0.997 | 0.997 | 0.994 | 0.999 |

$D$ | 0.551 | 0.581 | 0.538 | 0.584 |

${\sigma}^{2}$ | 0.997 | 0.998 | 0.995 | 0.993 |

${\rho}_{N}$ | 0.555 | 0.560 | 0.516 | 0.589 |

All | 0.998 | 0.998 | 0.996 | 0.999 |

## Table 10

Performance of the system for the extracted features in terms of Az values for tissue data (mouse colon in vitro) for classification between case Mouse 2 versus Mouse 3.

500.33nm | 550.02nm | 575.02nm | 625.24nm | |
---|---|---|---|---|

${N}_{c}$ | 0.541 | 0.596 | 0.596 | 0.571 |

$E$ | 0.938 | 0.961 | 0.974 | 0.906 |

$D$ | 0.556 | 0.578 | 0.509 | 0.622 |

${\sigma}^{2}$ | 0.816 | 0.896 | 0.887 | 0.784 |

${\rho}_{N}$ | 0.573 | 0.591 | 0.501 | 0.615 |

All | 0.940 | 0.968 | 0.987 | 0.913 |

## 5.

## Conclusions

A textural model that models the reflected scattered signal is presented, and its sensitivity to simulated changes in the morphology and physical characteristics of mimicking tissue being optically scanned is presented as a proof-of-concept to model scattered signals reflected from mucosal tissues. At the core of that method is the modeling of the scattered light using a SDM texture model that captures how the tissue consisting of diffuse and specular scatterers reflects light shined at different wavelengths to yield textured data and corresponding image model parameters (obtained through the SDM) that correlate closely with these changes in the tissue morphology and physical characteristics (for example in dysplasia formation). In addition, the method statistically analyzes 2-D scans of light scattering data, hence allowing for the assessment of correlation and textural characteristics that otherwise could not be revealed when the analysis of the scattering signal is a function of wavelength or angle, which has been done by other investigators. This proof-of-concept is shown using simulation, phantom data, and on a limited preliminary *in vitro* animal experiment to track mucosal tissue inflammation over time. Although the scope of the animal experiment is limited, it tests the applicability of the model on real tissue. The test scenarios for the simulation and the phantom examples are driven by biological applications such as dysplasia formation and detection. The tissue experiment is geared more toward inflammation tracking, and the ability of the SDM to differentiate between the various stages of inflammation as it develops in the mouse colon epithelium. This is of particular relevance, as there is a functional relationship between inflammation and cancer. The performance is measured by computing the area
${A}_{z}$
under the ROC curve for each parameter separately, and by fusing all the estimated parameter sets together. Combination of all the features results in an
${A}_{z}$
value up to 1
$({A}_{z}>0.991)$
(perfect detection for any false alarm rate) for simulated data, and
${A}_{z}>0.927$
for the phantom data. The performances for differentiation between pairs of various levels of inflammation are 0.859, 0.983, and 0.999. These encouraging results show that the proposed texture model and technique are excellent candidates to be applied to biological tissues, with the potential for enhanced endoscopy for detection of changes in epithelial tissue.

## Acknowledgments

We would like to thank Drexel Photonics Laboratory members Elina Vitol, Timothy Kurzweg, and Bahram Nabet for providing the optical device (probe, spectrometer, and light source) and the phantoms used in this research. Special thanks are also due to Elina Vitol for help with the phantom data collection. We also would like to thank Amol Karwa for his help during tissue experiments and for sharing useful biological information.

## Appendices

### Appendix: Mie Scattering and Signal Formation for Stochastic Simulator

The Mie scattering theory describes the light scattered by spherical particles.^{57} Scattering of light in tissue mainly consists of two components: singly scattered and multiply scattered. Singly scattered light, as predicted by the Mie theory, is not randomized and contains information about individual scatterers. On the other hand, multiply scattered light is thoroughly randomized. Diffusely scattered light from tissue contains information about its basic structures. In the epithelial tissue, the light transport is dominated by multiple scattering.^{58, 59, 60}

For multisphere scattering simulation, we used the T-matrix approach. The main complication here is how the presence of another scatter center in the neighborhood affects the radiation patterns, scattering matrix, and efficiencies. When two or more identical spheres aggregate into a cluster, the resulting composite particle is nonspherical.^{61} The Mie theory, although widely used in the literature as a scattering approximation theory, requires, modeling a cell as a homogeneous sphere.^{62} Various techniques have been reported for generalizing the Mie theory^{62} to deal with the application of it to a volume of spheres with various size distributions,^{63} and for modeling a spherical cell containing a homogeneous spherical nucleus.^{64} Some used anomalous diffraction approximations,^{65} some multiple solutions,^{66} some used 3-D finite-difference time-domain (FDTD) models of cellular scattering,^{45} and finally some used T-matrix computations.^{67} The problem of dependent scattering^{68} in the presence of a cluster of spheres has been extensively investigated in the literature.^{69, 70, 71, 72, 73} All of these reported techniques have important advantages over the Mie theory and helped in alleviating some of its limitations.^{62} Recent work has verified the practical applicability of the T-matrix method to clusters of spheres.^{74} The T-matrix method surpasses other frequently used techniques in terms of its ability to deal with a range of sphere sizes, as well as to systematically deal with nonspherical scattering based on the calculation of thousands of particles in random orientation.^{75} It is for these reasons that we chose the T-matrix approach in our multisphere scattering simulation experiment.

The strategy of the T-matrix solution of the direct light scattering problem can be introduced with the following algorithm: expanding the incident field in vector spherical wave functions regular at the origin; expanding the scattered field outside a circumscribing sphere of scatterers in vector spherical wave functions regular at infinity; and transforming the expansion coefficients of the incident field into those of the scattered field via a T-matrix.^{76} The relationship between incident and scattered electric field components perpendicular and parallel to the scattering plane as observed in the farfield is described by the amplitude scattering matrix:

## Eq. 10

$$\left[\begin{array}{c}{E}_{\simeq s}\\ {E}_{\perp s}\end{array}\right]=\mathrm{exp}\phantom{\rule{0.2em}{0ex}}\frac{\left[i\mathbf{k}(r-z)\right]}{-i\mathbf{k}r}\left[\begin{array}{cc}{S}_{2}& {S}_{3}\\ {S}_{4}& {S}_{1}\end{array}\right]\left[\begin{array}{c}{E}_{\simeq i}\\ {E}_{\perp i}\end{array}\right],$$^{80}that specifies the set of parameters associated with the phase and polarization of radiation. For a cluster of spheres, the relationship between incident and scattered Stokes parameters is described by the Mueller scattering matrix

^{57, 61}as:

## Eq. 11

$$\left[\begin{array}{c}{I}_{s}\\ {Q}_{s}\\ {U}_{s}\\ {V}_{s}\end{array}\right]=\frac{1}{{k}^{2}{r}^{2}}\left[\begin{array}{cccc}{S}_{11}& {S}_{12}& {S}_{13}& {S}_{14}\\ {S}_{21}& {S}_{22}& {S}_{23}& {S}_{24}\\ {S}_{31}& {S}_{32}& {S}_{33}& {S}_{34}\\ {S}_{41}& {S}_{42}& {S}_{43}& {S}_{44}\end{array}\right]\left[\begin{array}{c}{I}_{i}\\ {Q}_{i}\\ {U}_{i}\\ {V}_{i}\end{array}\right],$$^{77, 78, 79}) such as:

^{61}