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\text{to}+20\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\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\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\mathrm{deg}$ ( $\theta$ value) and $90\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\mathrm{deg}$ and $\phi =90\mathrm{deg}$ [Eqs. 12, 13]. The intensity is calculated from the scattering matrix using Eq. 14.

Fig. 2

Distribution of scattering spheres within 256 uniformly distributed cluster ensembles with each cluster having $N$ spheres each with diameter $d$ .

To introduce randomness into Mie scattering,¹⁶ 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:

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\mu \mathrm{m}$ diameter. For the simulation of the diseased cells, each cluster ensemble is considered to consist of four cells of $12\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\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\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.

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\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\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\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.

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.⁵⁴ 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.⁵⁴ 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\mathrm{nm}$ . We examined microspheres at diameter sizes of 3 and $10\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\mu \mathrm{m}$ .⁵⁴ The probe has a special $30\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\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.

Fig. 3

3-D schematic figure of the experimental setup used for data collection. The probe is held and moved by the motorized micromanipulator ( $xyz$ stage) allowing 2-D scanning of the sample. It is connected to a light source and spectrometer.

Fig. 4

Real pictures of the experimental setup used in data collection. This corresponds to the schematic figure in Fig. 3.

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\mu \mathrm{m}$ . Three 1-D scans are taken for each phantom. The scans are taken at 0, 2250, and $3500\mu \mathrm{m}$ (the $y$ direction). For each 1-D scan, 51 points are taken $90\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.⁵⁵ 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.⁵⁶ 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\text{to}8\min$ to complete the scan. A schematic of the animal study is shown in Fig. 5 .

Fig. 5

Animal study. Colon was excised on days 4 and 7 of 4% dextrin sulfate solution (DSS) feeding, causing increased inflammation with feeding. The excised colon was mounted in a petri dish containing phosphate buffered saline and was imaged within thirty minutes after mounting.

Three 1-D scans of the entire piece of the colon $(\mathrm{10,000}\times 3000\mu \mathrm{m})$ were taken in a zigzag fashion, using a micrometer stage within $30\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\mu \mathrm{m}$ . Each scan consists of 50 points. The step size is taken as $200\mu \mathrm{m}$ . After the data are collected, the same processing steps as with the phantom data are applied.

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.

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.

Fig. 6

3-D plots of extracted features in 3-D feature space. The features are: $D$ (Rayleigh scattering degree of the diffuse component), ${\sigma }^2$ (residual error variance of the diffuse component) and ${\rho }_N$ (normalized correlation coefficient of the diffuse component). These are shown for the simulation pairs: (a) L1-S1, (b) L2-S2, (c) L1-L2, and (d) S1-S2. (refer to Table 1). The good separation between different simulated sets can be observed from the given plots.

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.

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.

Fig. 7

3-D plots of extracted features in 3-D feature space. The features are: $D$ (Rayleigh scattering degree of the diffuse component) ${\sigma }^2$ (residual error variance of the diffuse component), and ${\rho }_N$ (normalized correlation coefficient of the diffuse component). These are shown for the simulation pairs: (a) A1-A2, (b) A2-A3, and (c) and A1-A3 (refer to Table 2). The good separation between different simulated sets can be observed from the plots.

Table 5

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

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\text{to}600\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.

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.

4.3.

Performance Evaluation on Tissue

The classification is reported for wavelengths between $450\text{to}650\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\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\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\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\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.

Fig. 8

Scatter plot of mean energy of coherent scatterers versus residual error variance for the three animals studied (Mouse1, Mouse2, and Mouse3). The scatter plot shows good separation of extracted features between the three levels of inflammation. As the disease level increased, so did the values of these two features.

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.

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.

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.