Differing self-similarity in light scattering spectra: A potential tool for pre-cancer detection

The fluctuations in the elastic light scattering spectra of normal and dysplastic human cervical tissues analyzed through wavelet transform based techniques reveal clear signatures of self-similar behavior in the spectral fluctuations. Significant differences in the power law behavior ascertained through the scaling exponent was observed in these tissues. The strong dependence of the elastic light scattering on the size distribution of the scatterers manifests in the angular variation of the scaling exponent. Interestingly, the spectral fluctuations in both these tissues showed multi-fractality (non-stationarity in fluctuations), the degree of multi-fractality being marginally higher in the case of dysplastic tissues. These findings using the multi-resolution analysis capability of the discrete wavelet transform can contribute to the recent surge in the exploration for non-invasive optical tools for pre-cancer detection.


Introduction
The use of optical techniques for the study of biomedical systems is a rapidly developing field that has seen a dramatic expansion in the recent years, partly due to tremendous progress in the field of lasers, fiber optics and associated technologies. Both medicine and biotechnology require appropriate instrumentation to analyze and monitor biological systems for deviations from normality. Optical methods, due to their non-invasive nature, are providing novel approaches for medical imaging, diagnosis and therapy. Considerable efforts have been made in the recent past to use optical techniques such as fluorescence spectroscopy [1,2,3,4], Raman spectroscopy [5] and elastic scattering spectroscopy [6] for quantitative and early diagnosis of various diseases. Several optical imaging techniques like coherence gated imaging, polarization gated imaging and diffuse optical tomography are also being actively pursued for obtaining high resolution (micron scale) images of biological objects and their underlying structure [7,8,9,10,11].
For optical diagnosis, elastic and inelastic light scattering spectra from tissues are exploited. The in-elastically scattered light (via processes like fluorescence and Raman) contain useful biochemical information about the sample that can be employed for probing subtle biochemical changes as signatures of disease progression. On the other hand, elastically scattered light from biological tissues contain rich morphological and functional information of potential biomedical importance [6,12,13,14,15,16,17,18,19,20,21,22,23,24,25]. Both the angular and wavelength dependence of the elastically scattered light from tissue can be analyzed to extract and quantify subtle morphological changes taking place during progression of a disease [14,15,16,17,18,21,22], and thus may be explored as a sensitive tool for early diagnosis. This would however involve appropriate modeling of light scattering in complex random media like tissues, and the development of suitable approaches to extract/interpret the morphological information contained in the elastic light scattering signal.
The spatial fluctuation of the refractive index in biological tissues arising from scatterers ranging in sizes from a few nanometers to a few micrometers give rise to elastic scattering [6,20]. The lack of sufficient knowledge about the complex dielectric fluctuations in the tissues pose a formidable problem in the exact modeling of light scattering. Nevertheless, several efforts have been made in the recent years using electromagnetic (EM) theory based approaches like Mie theory and Born approximation to model and understand the scattering process from biological tissues [26,27,28,29,30]. It has also been shown that the refractive index fluctuations in biological tissues are fractal in nature which can be used to understand the structural changes in tissues induced by diseases [16,27,31,32,33].
Since the tissue morphology dependent refractive index fluctuations are recorded in the elastic scattering spectra [34], analysis of elastically scattered spectral fluctuations using sophisticated fluctuation analysis tools might facilitate extraction and quantification of subtle morphological changes associated with early stages of cancer. The scaling behavior which is generally assumed to be global (mono-fractal), has been shown to manifest in the local fluctuations in various physical processes [35,36] and has been characterized using Multi-Fractal De-trended Fluctuation Analysis(for example see [37]). Wavelet Based Multi-Fractal Detrended Fluctuation Analysis (WB-MFDFA) is one other state-of-the-art technique that can be used for extracting and quantifying the self similarity at varying length scales associated with the structural changes associated with cancer progression due to the inherent use of fractal like transformation kernels.
Wavelet transform due to it's multi-resolution analysis capability using the Daubechies' basis which extract the polynomial trends (for example, Db-4 and Db-6 extract the linear and quadratic trends respectively) has been shown to characterize the scaling behavior and selfsimilarity of empirical data sets quite faithfully [38,39]. Indeed, it has been initially explored to analyze tissue fluorescence spectra in an attempt to distinguish between normal and dysplastic tissue [40,41,42,43,44]. In this work, we employ this multi-resolution property of wavelets to ascertain the changes in the self-similarity of dysplastic human cervical tissues as opposed to healthy human cervical tissues by analyzing the esoteric nature of the fluctuations in tissue light scattering spectra.
This article is organized as follows: Fourier and power spectrum analysis is reviewed in 2.1, discrete wavelet transform in 2.2, wavelet based power law analysis in 2.3, wavelet based multi-fractal de-trended fluctuation analysis in 2.4 and correlation based analysis in 2.5. Sec. 3 describes the experimental methods for light scattering measurements from tissues. Sec. 4 deals with our findings from the analysis and contains a discussion of the same in the context of the differences between the normal and dysplastic samples. In Sec. 5 we conclude with the prospect of pre-cancer detection using light scattering techniques combined with novel fluctuation analysis methods.

Fourier analysis and power law spectrum
Fourier Analysis has traditionally been a preferred tool for analysis of experimental data sets. Here, we just briefly review the Discrete Fourier Transform (DFT) and its power spectrum. For a data set x(n), n = {1, 2, . . .}, the DFT is a linear transformation over an orthogonal basis given by: and it's power spectrum is known to follow a power law behavior for self-similar processes The power law coefficient α is related to the Hurst scaling exponent H ∈ (0, 1), and the fractal dimension D f which identify the nature of fractal behavior by α = 2H Though this analysis assumes a global scaling behavior, due to the presence of inhomogeneities in the intra-cellular structure, the scaling behavior turns out to be localized in the spatial frequency domain. In multi-fractal analysis, a local Hurst exponent is calculated which quantifies the local singular behavior and provides useful and minute information which is hidden between different scales. We therefore applied a more general wavelet based approach for the multi-fractal analysis.

A brief review of wavelet transforms
Wavelet Transforms [45,46,47,48], in both the forms of Discrete and Continuous transforms in the recent years have emerged as an invaluable tool in the field of data analysis and interpretation. Here, we will briefly review the Discrete Wavelet Transform (DWT).
In DWT, two functions, namely, φ (n) and ψ(n), called the father and the mother wavelets respectively, form the kernels. They satisfy the admissibility conditions: The scaled and translated versions of the mother wavelet ψ(n) are called the daughter wavelets which form a complete set and differ from the former in terms of their height and width. At the s th scale, the height and width of the daughter wavelet are 2 s and 2 s/2 of that of the mother wavelet respectively. s and m are the scaling and translation parameters. The DWT for a function f (t) is given by, the coefficient a m (called the approximation coefficient) extracts the trend and d s,m (called the detail coefficient) extracts the details or fluctuations from the signal. Here, l = ⌊log(N)/ log(2)⌋ is the upper bound for taking the maximum scale for analysis above which the edge effects corrupt the wavelet coefficients giving rise to spurious results. This mathematical artifact is explained by the cone of influence [47].

Wavelet Fluctuation Power Law Analysis
The Daubechies' family of wavelets are made to satisfy vanishing moments conditions, and hence are able to isolate polynomial trends from a time series (see [45]). The different wavelets in this family isolate trends of different polynomial orders, for example, Db − 4 isolates a monomial trend while Db − 6 isolates a quadratic trend. We first obtain the profile from the fluctuations by taking their cumulative sum: where, Y and X are profile and fluctuation signal respectively and N is the data length. In this method, we first obtain the fluctuations at every scale by a wavelet reconstruction taking the approximation coefficients (using Db-4 wavelet) and subsequent subtraction from the signal. A flowchart for this fluctuation extraction is shown in Fig. 1. We then fit the Fourier power spectrum of these fluctuations to obtain the power law exponent as a function of scale, i.e. α ≡ α(s) for the identification of the short and long term correlations. This technique sheds light on the scaling behavior of the fluctuations in different spatial frequency regimes.

Wavelet Based Multi Fractal De-trended Fluctuation Analysis
The Multi Fractal De-trended Fluctuation Analysis (WB-MFDFA) algorithm proposed by Manimaran et al. [38,39] has been employed gainfully to extract the multi-fractal nature of a variety of physical systems. In order to apply the WB-MFDFA algorithm, we then obtain the fluctuations using the extraction algorithm shown in Fig. 1. The asymmetric nature of the wavelet function and the edge-effects (due to the cone of influence) encountered during the convolution can affect the precision of the fluctuations. Hence, we repeat this procedure on a reversed profile and then take the average to get the fluctuations at every scale which are denoted by Fs. We useds to represent the wavelet scale so as not to confuse with the segment length s which is related to the wavelet scales by the number of filter coefficients for a given wavelet. These fluctuations are then segmented into M s non-overlapping sections such that M s = ⌊N/s⌋, where, s and N are the window size and the length of the fluctuations respectively. Subsequently, we obtain the q th order fluctuation function F q (s) by where the order of moments q can and both positive and negative integers. The fluctuation function for self similar processes follows a scaling law, the scaling function given by F q (s) ≈ s h(q) . We should note here that the smaller fluctuations in the light scattering spectrum will be influenced by the negative q values, whereas the positive values of q will impact the larger fluctuations. The scaling function h(q) calculated at q = 2 corresponds to the Hurst exponent [35]. H = 0.5 represents uncorrelated (white noise, f 0 ) or brown noise f −2 , while H > 0.5 represents long range correlations or persistent behavior. H < 0.5 reveals short range correlations or antipersistent behavior. Mono-fractals are scale independent and hence their h(q) is independent of q. They can be characterized by a single parameter like the fractal dimension. However, for multi-fractals, the h(q) is not independent of q and they require a more complex function like the singularity spectrum for its characterization [36]. The classical multi-fractal scaling exponent τ(q) defined by standard partition function based formalism [50,51] is related to the Hurst exponent h(q) by τ(q) = qh(q) − 1. The multifractality can also be characterized by using the singularity spectrum f (β ), which is related to τ(q) through a Legendre transform: β = d/dq[τ(q)] and f (β ) = qβ − τ(q); where β is the singularity strength or Hölder Exponent [37]. f (β ) denotes the dimension of the subset of the series to be characterized. β , f (β ) and h(q) are related by β will have a constant value for mono-fractal series, whereas for multi-fractal series β values will have a distribution. Just like in the case of h(q), a constant τ(q) will signify a monofractal, while a multi-fractal will depend on the order of the moments. The singularity spectrum is broader when correlations are present in the series and if the correlations are absent or lost, then the singularity spectrum becomes narrower.

Correlation based analysis
The correlation matrix analysis method has been extensively used to study various physical and biological systems (for example see [52]). For two variables x i ∈ X and y j ∈ Y , the correlation matrix is given by: where, . . . is the mean. It is well known that multi-fractality is related to the domain structure formation in the correlation matrices and we present the correlation matrix analysis in the wavelength domain to support the wavelet based multi-fractal de-trended fluctuation analysis.

Experimental materials and methods
Pathologically graded (CIN-I or dysplastic) unstained sections of human cervical tissue and their normal counterparts were frozen and vertically sectioned (thickness ∼ 5µm, lateral dimension ∼ 4mm × 6mm) on glass slides were used for light scattering measurement and analysis. A schematic representation of the experimental setup is shown in Fig. 2. The angle resolved elastic scattering spectra were recorded using a goniometric arrangement for the spectral range λ : 400nm − 800nm and angular range θ : 10 • − 150 • at 10 • intervals. White light output from a Xe-lamp (Newport USA, 50-5000 W) was collimated using a combination of lenses and were made incident on the sample kept at the center of the goniometer. The spot-size incident on the sample was controlled by a variable aperture (∼ 1mm). The scattered light from the sample was collimated by a pair of lenses and was then focused into a collecting fiber probe, the distal end of which was coupled to a spectrometer (Ocean Optics HR2000). The recorded spectra at each scattering angle was normalized by the lamp spectra.

Results and Discussions
Typical elastic light scattering spectra recorded from normal and dysplastic tissues are shown in Fig.3. Representative spectra for forward and back-scattering are shown for θ = 10 • and The sample mounted on a goniometer is illuminated by a white light collimated from a Xe-source and the scattered light is then collimated using two lenses (L 2 and L 3 ) and recorded using a Fiber-Optic Spectrometer. The θ f represents the forward scattering angle, while θ b is the backward scattering angle. The forward and backward scattering process are represented by red and blue lines respectively while the Fiber-Optic Spectrometer setup is shown in green. 0 ≤ θ ≤ 180 • is the scattering angle at which the spectra were recorded. We must note in the gray shaded region extending from 80 • to 110 • , the scattering spectra was not recorded due to geometry induced artifacts. . We can see that the spectrum flattens as we move towards the backward scattering angles as compared to the forward scattering angles, which signifies the involvement of larger sized scatterers, such as the nucleus, in forward scattering and smaller sized intracellular organelles in the backward scattering. Such angular dependence of elastic scattering spectra from tissues have been extensively studied both experimentally and theoretically. These studies have indeed shown that the contribution of larger scatterers (like cells, nuclei) are more dominant in the spectra recorded at forward angles, whereas the spectra recorded at back-scattering angles are typically more influenced by the smaller sized scatterers [24]. Although it is difficult to make a one-to-one correspondence between the intensity fluctuations and local refractive index fluctuations in the morphological structures, some information about the morphology can still be inferred from the light scattering spectral fluctuations. Indeed, it has been shown that with Born approximation, the light scattering spectrum can be represented as a Fourier transform power spectrum of the local refractive index fluctuations (see Sec. 2.1 for details). Thus, the observed spectral fluctuations contain information about the corresponding refractive index fluctuations in the Fourier domain as discussed below.
The results of the Fourier analysis on the light scattering spectra for different scattering angles is shown in Fig. 4. In the forward scattering region (10 • − 50 • ), we observe that the dysplastic samples show ∼ 0.7 ≤ α ≤ 1.2, in the same region, however, the normal tissues show 0.4 ≤ α ≤ 0.8 (where α is the power law coefficient). The higher values of the powerlaw coefficient α corresponds to higher values of the Hurst scaling exponent H, indicative of   "coarseness" of the fluctuations tending towards sub-fractal behavior; whereas at higher angles (110 • − 150 • ), α values are observed to be smaller in dysplastic samples than that for normal samples in the same region indicating more "roughness" in the fluctuations, signifying a trend towards extreme fractality.
The higher α values in the angular range (10 • ≤ θ ≤ 70 • ) for the dysplastic tissues as shown in Fig. 4, due to the more dominant contribution of the large scale scattering structures. This possibly can be related to the proliferative nuclear morphology in dysplasia [13]. This follows because such morphological changes are expected to lead to coarser fluctuations in light scattering spectra as is evident from the higher value of H in the dysplastic tissues in the forward scattering range. On the other hand, the lower values of α in the back-scattering angles (110 • ≤ θ ≤ 150 • ) indicates the predominance of small scale heterogeneity in the dysplastic tissues. It is also pertinent to note that around 60 • , there is no distinct difference, possibly due to the simultaneous contribution of different sized scatterers in this angular range resulting in "lumped" effects.  In Fig. 5, we present the plots of power-law scaling exponent α as a function of the scale for normal (Fig. 5(a)) and dysplasia (Fig. 5(b)) for forward scattering angles 40 • and 50 • and for backward scattering angles 140 • and 150 • . This analysis was performed following the discussion in Sec. 2.3. It must be noted that the difference in the absolute value of α for Fourier and wavelet analysis arises from the fact that, in Fourier analysis, we analyze the signal itself while in wavelet analysis, we obtain the α values from the fluctuations which can be thought of as higher derivatives of the signal. We observe that the α values show a strong dependence on the scattering angle θ . In both Fig. 5(a) and Fig. 5(b), we observe the presence of a broad spectrum of processes such that 1.4 ≤ α ≤ 3.2. This arises possibly from the contribution of differently sized scatterers as has been mentioned earlier. We must however note that the behavior of the normal and dysplastic tissues in power law regime turn out to be very different. For example, in the forward scattering region, at 50 • , we observe that the normal tissues show an average α = 2.6, while the dysplastic tissues show an average α = 2.0 over all scales. Similarly, in the backward scattering region, we obtain α = 1.8 and α = 3.1 for normal and dysplastic tissues respectively at 150 • . In the other regions also, we find stark differences between normal and dysplastic samples in the α values averaged over all scales. As discussed in Sec. 2.4, we have depicted the Hurst exponent as a function of the scattering angle θ in Fig. 6; where we notice that though for both the normal and dysplastic tissues show a "near-random" behavior corresponding to H = 0.5 (H = 0.5 is indicative of Brownian behavior), we observe that the dysplastic samples show a marked departure from the normal samples in their trend. As observed for α(θ ) in Fig. 4, H has a higher value in the forward scattering range 10 • ≤ θ ≤ 70 • while the trend reverses in the backward scattering range 120 • ≤ θ ≤ 150 • .
The singularity spectrum f (β ) is a quantitative indicator of the exact nature of the selfsimilarity and its width represents the strength of the multi-fractality. In Fig. 8, we have shown the f (β ) as obtained from Eq. 7 as discussed in Sec. 2.4. We observe that the dysplastic samples have a higher multi-fractality than the normal samples, indicated by the width of singularity spectrum. It must be noted that for a mono-fractal, the singularity spectrum is similar to a Gaussian with a very small variance. This is consistent with our observations of Fig. 7, where, the plots of h(q) vs q for normal and dysplastic tissues at a few representative forward and backward scattering angles are shown. For mono-fractals, the h(q) is independent of q. In Fig. 7, we observe that for the normal samples, at angles 50 • − 70 • , the dependence of h(q) on q is very small (∆h(q) ∼ 0.2), which implies a trend towards mono-fractality. The normal samples thus show a wide range of variation from mono-fractals to multi-fractals. However, for dysplasia,  the h(q) dependence on q is high (∆h(q) ∼ 0.7 ) for all scattering angles, which indicates that the dysplastic samples show more multi-fractality than the normal samples. The spectral correlation matrices are shown in Fig. 9 (following Sec. 2.5). It is clear from Fig. 9(a) and 9(b) that though the normal tissues do not show any distinct correlation sectors in this domain, other than the expected correlation that would occur at neighboring wavelengths; dysplastic tissues show the presence of three dominant sectors. It must also be noted that the range of correlation increases from 0.70 − 1.00 for normal to 0.97 − 1.00 for dysplasia. This indicates a higher correlated behavior of the dysplastic tissues than the normal tissues in addition to domain formations in the spectral range. This possibly arises due to the fact that during dysplastic progression, the homogeneous cell morphology gives way to a more fragmented and heterogeneous structure.

Conclusion
In conclusion, we have applied a combined Fourier based and discrete wavelet based analysis on the fluctuations extracted from the elastic light scattering spectra of normal and dysplastic human cervical tissues. This approach clearly revealed otherwise hidden signatures of self-similarity in spectral fluctuation for both normal and dysplastic tissues, with significant differences in the nature of self-similarity. Fourier analysis of these fluctuations indicated the existence of multi-fractal nature and was further confirmed by WB-MFDFA. Dysplastic tissues showed marginally higher multi-fractality over the entire angular region compared to their normal counter-parts. The scaling exponent was observed to have angular dependence, possibly arising from the size distribution of the scatterers present in such complex systems. Note that this novel fluctuation analysis approach was initially explored on elastic scattering spectra recorded from a limited number of tissue samples, and the results were qualitatively similar. A more systematic study on statistically significant number of samples is currently underway, and the results will be reported in the near future. Never-the-less, early indications show promise of this approach for quantification of the morphological alterations associated with pre-cancers.