Open Access
Add to BibSonomyAbstract
A novel Fourier-based image analysis method for measuring fractal features is presented which can significantly reduce artifacts due to non-fractal edge effects. The technique is broadly applicable to the quantitative characterization of internal morphology (texture) of image features with well-defined borders. In this study, we explore the capacity of this method for quantitative assessment of intracellular fractal morphology of mitochondrial networks in images of normal and diseased (precancerous) epithelial tissues. Using a combination of simulated fractal images and endogenous two-photon excited fluorescence (TPEF) microscopy, our method is shown to more accurately characterize the exponent of the high-frequency power spectral density (PSD) of these images in the presence of artifacts that arise due to cellular and nuclear borders.
©2012 Optical Society of America
1. Introduction
The analysis of biomedical images is critical for detection of abnormalities and disease, but it is often subject to the interpretation of a medical professional. Starting as early as the 1960s, efforts have been made to develop quantitative tools based on automated image analysis algorithms to assist physicians and researchers in characterizing tissue properties [1]. Optimization and development of these methods is still underway, and more groups are recognizing the utility of these techniques for extracting patterns and information from biomedical images. Uncovering this image information is likely to lead to the discovery of novel and objective diagnostic criteria, improving diagnostic sensitivity and enabling earlier disease detection. In this way, widespread application and optimization of quantitative image analysis techniques has great potential to impact the performance of clinical diagnostics and basic research that relies on interpretation of biomedical images.
Fourier-based techniques have wide-range applications in signal and image assessment and are gaining a more critical role in tissue characterization. For example, these techniques have been implemented to characterize bone structure in computed tomography images [2] and to detect cardiac arrhythmias in electrocardiogram signals [3]. The squared amplitude of the Fourier transform (FT) is referred to as the power spectral density (PSD). Biomedical images of cells and tissues often exhibit PSDs with inverse power-law frequency dependence (i.e., proportional to k-β, where k is spatial frequency and β is the power-law exponent), which can indicate a scale-invariant (fractal) organization of the imaged features [4]. Scale-invariance describes features or patterns that persist over multiple length scales. These features must satisfy conditions of self-similarity; meaning a fractal object is similar to a subset of itself. In special cases, fractal can be considered self-affine if variation in one direction scales differently than variation in another direction [4–6].
Fractals are present in a wide variety of natural systems [5, 6], including rock strain distributions [7], microvascular networks [8], and chromatin aggregation [9]. Numerous groups have reported on the fractal nature of subcellular inhomogeneities and their variation with disease state [9–16]. However, many of these studies are based on light scattering properties of cells and tissues, providing an indirect (and non-singular) determination of cell and tissue morphology. Quantitative characterization of fractal features can thus vary with the particular model assumed for subcellular or tissue density fluctuations; for example, whether these features exhibit von Kármán [11, 13], exponential [10, 17] or stretched exponential spatial correlations [16]. Furthermore, direct analysis of subcellular features is often attained by invasive methods, such as by histological staining [17] or electron microscopy [18], which may alter morphology from its nascent state.
Biomedical images of living cells or tissues obtained via fluorescence microscopy, on the other hand, have the advantage of providing a direct and noninvasive means for determining cell morphology down to submicron length scales. In this study, we characterize images that rely on intrinsic fluorescence from nicotinamide adenine dinucleotide (NADH), which emanates from cell mitochondria [19]. Mitochondria are the main energy-converting organelles of mammalian cells [20]. Metrics that quantitatively assess the degree of correlation of mitochondrial networks could serve as useful indicators of cellular health status. For example, early work by Hackenbrock showed that mitochondrial networks rapidly become more condensed in liver cells in which oxidative phosphorylation was activated [21–23]. More recent work has focused on the thinning and branching of the mitochondrial networks upon the switch from glycolytic to oxidative energy metabolism, which is thought to be a critical indicator of cell growth or differentiation [24, 25].
For PSD-based analysis, as the measured power-law exponent (β) of an image increases, the fractal metrics we assess indicate the mitochondrial networks become more correlated. Generally, when the PSD spectral content is flat (β = 0), the signal is uncorrelated, white Gaussian noise (Table 1 ). Images that contain features characterized by power exponent values ranging between 0 < β < 2 are considered fractional Gaussian noises (fGNs), whereas fractional Brownian motions (fBMs) are characterized by exponent values between 2 < β < 4 [4]. fBMs are distinct from fGNs in that they have long-range correlation (“memory”), which has been described as a statistical relation between the increments of data set values [4]. Visually, these long-range correlations result in progressively more clustering as the power-law exponent increases.
Table 1. β values and corresponding statistical processes
In practice, most PSD-based analyses that assess the fractal nature of biomedical images do not account for artifacts that can be produced by edge effects from cellular and nuclear borders [14, 26]. To overcome these artifacts, some studies have focused on small regions inside cells, which requires sub-image selection and drastically limits resolution and field of view. Others have discussed in a geological context how power spectral estimates of scale-invariant data sets depend on factors other than the fractal characteristics of the given data set, such as sampling, noise, and, edge effects [27].
The objective of our study is to develop an automated, PSD-based analysis for the quantitative characterization of internal morphology (texture) of image features with well-defined borders that overcomes the limitations introduced by the presence of these borders, in order to develop more robust diagnostic biomarkers of organizational attributes. We present a mechanistic computational study in which we evaluate simulated images of scale-invariant (fractal) patterns combined with controlled cell-shaped features to show that these cell-shaped features and image background have predictable effects on PSD power-law decays. Based on this evaluation, we quantify the manner in which the presence of cell- and nuclear-shaped features hinders the ability of PSD-based methods to accurately assess intracellular patterns. To overcome this limitation, we have developed a digital object cloning (DOC) method for use prior to traditional PSD analysis. We characterize simulated images that contain structures and textures common to epithelial tissues and compare the ability of PSD-based methods to accurately quantify intracellular texture both with and without DOC pre-processing. Finally, to showcase the diagnostic utility of this improved technique, we apply it to characterize experimentally acquired autofluorescence images from engineered normal and pre-cancerous epithelial tissues
2. Methods
2.1 Simulations of scale-invariant images
To test the ability of PSD-based approaches to sense changes in fractal scaling, simulated fractal images, , of size N × N and which have power spectral density that scales as a power-law with pixel values represented by the vector , , are generated in MATLAB based on Voss’s inverse Fourier filtering method [28]. In order to construct a fractal image using this approach, white noise, , is generated with a mean of 0 and a standard deviation of 1 with the randn function in MATLAB. The discrete Fourier transform, , as a function of radial spatial frequency, , where and for is taken as:
Next, the phase, , of the Fourier transform of is computed by dividing the Fourier transform (real and imaginary parts) by its magnitude. For :
The phase is filtered with a transfer function, , to produce the random function,
which is retransformed into the spatial domain to produce a fractal image, , where
To generate inverse power-law dependent PSD scaling from requires that:
where β is the exponent of the power-law decay and for β values ranging from 0 to 4 [28]. For each generated image, Eq. (6) is fit to the radial (angularly-averaged) power spectral density (PSD):
where is the fit to the radial PSD and is the radial spatial frequency. To allow direct comparison to acquired two-photon excited fluorescence (TPEF) images, radial spatial frequency of simulated images is reported in inverse microns (two pixels correspond to one micron in length). A short, high frequency tail appears in all power spectra of images containing cell-shaped objects. To avoid fitting this tail, the maximum and minimum PSD values are determined. The high frequency region that contains values greater than the second percentile of the PSD amplitude are excluded from the fit. This high frequency limit (LHF), is used to restrict the fits to the frequency range 0.1 μm−1 < < LHF. This is the region in which inverse power-law behavior of the PSD spectra has been typically observed in previous work with autofluorescent images of epithelial tissues [14].
2.2 Engineered tissue constructs and TPEF data acquisition
Organotypic rafts cultured with normal human foreskin keratinocytes (HFKs) allow for generation of multilayered tissues that mimic the gradient of differentiation of native epithelium and are grown as described in detail previously [29]. Briefly, keratinocytes are plated on top of a neutralized bovine type I collagen (~4mg/ml) matrix with dermal fibroblasts, and raised to an air-liquid interface. The tissues are provided nutrients from below to simulate natural nutrient delivery by a vasculature bed for 10 days, at which point they are imaged. For pre-cancerous tissues, human-papillomavirus (HPV-) immortalized epithelial cells are used in place of normal HFKs as detailed previously [14]. TPEF images are acquired on a Leica TCS SP2 confocal microscope (Wetzlar, Germany) equipped with a Ti:sapphire laser (Spectra Physics, Mountain View, CA). Samples are placed on glass coverslips, excited with 755nm (TPEF) light and imaged using a 63x/1.2 NA water immersion objective, which yields 8-bit, 512 x 512 pixel images of 238 x 238 μm2, in approximately 1 s. TPEF images are acquired by a non-descanned PMT with a filter cube containing a 700 nm short pass filter (Chroma SPC700bp) a dichroic mirror (Chroma 495dcxr), and an emitter bandpass filter centered at 460 nm (Chroma 460bp40). This excitation and emission filter combination allows collection of fluorescence signal emanating primarily from mitochondrial NADH [14, 30, 31]. Five regions from three different tissue constructs for the diseased and healthy tissues are assessed with PSD analysis.
2.3 Creation of simulated cell objects (SCOs)
Images of simulated cell objects (SCOs) were created in order to characterize the effect of nuclear and cell border shape on the PSD power-law decay. To achieve this, binary images of randomly positioned circles on the size order of cells and nuclei are created in MATLAB that model spatial features of autofluorescence images acquired previously [32] and in this study. First, the MATLAB randi function is used to generate a 512x512 matrix with uncorrelated random values. To randomize circle location and control circle density, one value, which is repeated randomly throughout the matrix, is designated for circle location. To simulate nuclei, the same pixel locations are mapped and a second circle of smaller radius is drawn at the specified location. To investigate PSD sensitivity to a nuclear feature, 20 individual images are generated that contain 60, 35-micron diameter circles. This is repeated in a second set of simulations, which include a second smaller circle, the ‘nucleus,’ with a diameter of 12 microns. In biomedical images, variance in cell size is common, so a Gaussian distribution of circles with diameter variance equivalent to 2 microns is evaluated in every case, which is relevant to previous investigation of biomedical images of cells in a 3D collagen matrix [32].
The Fourier transform of a uniform circular object yields a Jinc function, which is related to a Bessel function of the first kind. According to the mathematical properties of the Jinc function, the distance from the first point to the first minimum in the Fourier spectrum will be 1.22/d, where the diameter of the circle is d [33]. To assess PSD sensitivity to circle size in our FT analyses, first minima of angularly-averaged PSD spectra, are found by locating the first zero crossing of the second derivative of the PSD.
In the high frequency limit, the PSD of a Jinc function is known to have oscillations that decay as an inverse power-law with a power exponent equal to 3. In order to examine the possible interference of the behavior of a Jinc function with the measurement of intracellular texture from cellular images, fractal images with known power-law decay values are digitally applied to the circular regions of the simulated images. To explore the sensitivity of this analysis to cell morphology and nuclear structures, we fit the PSD of the simulated images with Eq. (6) and the exponent of the fit is compared with input exponents of the fractal that has been digitally applied. To weaken contributions from background signal and edge effects, background is set to the average value of the intracircular signal when indicated in the text.
2.4 Identification of cell and nuclear borders for cloning TPEF images
In order to isolate the effect of cell and nuclear borders on PSD-based intracellular assessment, it is necessary to create a mask that separates intracellular regions from image background. Five TPEF images from each of the cell layers (superficial, para-basal, and basal) of three independent engineered epithelia from healthy and diseased groups are segmented so that only the cellular regions remain. Specifically, for each image, a set of binary images for every whole threshold value (0-255) for an 8-bit scale is produced in MATLAB. The power spectral density is then computed for each of the binarized images [34]. The binarization threshold value that creates an image with the highest peak power spectral intensity is defined as the optimal threshold to segment the cells from the background, because high peak power spectral intensity is related to the presence of low frequency image features. Finally, groups of pixels that consist of less than 80 interconnected pixels are removed. This process effectively creates a mask that isolates large objects, such as cells.
2.5 Creation of simulated fractal images with defined borders
Using the mask determined as described above, fractal images with known PSD power-law decay values (β = 0-4) are digitally applied to the intracellular regions according to this mask. The digital application of a fractal image is repeated four times, yielding a total of 60 simulated images with known intracellular texture (we start with five images from each of the three main epithelial tissue layers). To assess the contributions from the image background, a black background or the average value of the intracellular foreground is digitally applied to background regions. Power-law exponents are determined from fits to the power-law region of the radial PSD for all synthesized images. Measured exponents are then compared to input exponents from fractals that are simulated and digitally applied to the intracellular regions.
2.6 Digital object cloning (DOC)
We develop a randomized digital object cloning (DOC) method that isolates intracellular regions and clones these regions into the image background. This method automatically reduces artifacts related to image background as well as larger, periodic non-fractal features, such as cell size and intercellular tissue organization. Regions of signal are isolated using the aforementioned thresholding technique. In order to randomly fill the background of each thresholded image, two vectors with random values ranging from −512 to 512 are generated using the MATLAB randi function, which creates a matrix of uniformly distributed pseudorandom integers. Using the MATLAB function circshift, which circularly shifts the values in an array, a parallel image is created with the original foreground values shifted according to the values in the two vectors. The regions in the shifted, parallel image that overlap with background regions from the original image are added to the original image, filling the background and leaving the original signal in the foreground. This process is iteratively repeated until no original background pixels remain. Generally, five to ten iterations are needed to fill the image background.
2.7 Statistical testing and β error calculation
The correlation between measured and input β values for TPEF images is assessed via the Tukey Honestly Significant Different (HSD) test in JMP statistical software. Average β-error is calculated as the average of the absolute value of the difference between the input β value and the measured β value for each simulated image for β’s between 0 and 2.
3. Results
3.1 Simulated cell objects (SCOs) hinder the ability of PSD-based methods to accurately assess intracellular patterns
To demonstrate the baseline sensitivity of the PSD-based approach for recovering power-law exponents (β) that describe the fractal character of images, we evaluate the PSD of square fractal images generated as described in Methods. Figure 1 displays angularly-averaged radial PSDs (in log-log scale, Panel (a)) that correspond to simulated fractal images (Panels (b-f)), which become more clustered in gross-appearance with increasing power exponents (β). Measured power exponents, determined by fitting Eq. (6), match input exponents used in fractal simulations with negligible associated error (R2 > 0.99).
Fig. 1 (a) Power spectral density of generated scale invariant images of varying fractal character with fits (shown in panels b-f) for β = 0, 1, 2, 3, 4.
To assess whether variations in cell-shaped features impact PSD-based outcomes, we generated binary image models of simulated cell objects (SCOs) that consist of uniform circles on the size order of cells and nuclei. We digitally apply fractals generated with varying β values to model SCO images containing SCOs of 35μm diameter (Fig. 2(a) ), without nuclei, and a homogenous black background. Figure 2(b) displays the angularly-averaged radial PSDs for SCO images that are created with fractals that have β values of 1, 1.5, 2, 2.5, 3, and 4 as indicated. As input β value increases, the measured slope of the power-law decay (βM) trends toward the expected value. For values greater than 3, the characteristic decay of the Jinc function (βJinc = 3) dominates entirely over the fractal pattern decay. In all cases except β = 3, βM values are far from the true intracellular input β value, showing that the presence of SCOs hinders accurate intracellular texture characterization.
Fig. 2 (a) Split panel of two representative simulated cell object (SCO) images without (left) and with a nuclear structure (right) (b) PSD spectra from SCO images with digitally applied fractals and a black background showing there is no sensitivity to β values greater than 3. (c) PSD spectra from SCO images without (black) and with nuclear structures (gray) containing SCOs that have 35 μm and 12μm for cell and nuclear diameters, respectively, the average value of the foreground in the background, and a β = 2.8. A shoulder that can be attributed to the nuclear border can be seen in the spectrum of the model images with nuclear structures. The spectrum is smoother than in (b) as a result of the change in background value.
In biomedical images an inner nuclear structure is present within each cell and such images are generated as described in Methods (Fig. 2(c)). The addition of a nuclear structure affects the lower frequency region of the PSD spectrum as observed by the emerging shoulder in the representative spectrum (Fig. 2(c), gray arrow). The frequency region corresponding to this shoulder is related to the diameter of the nucleus, with a nucleus of 12 microns creating a local minimum at a frequency equal to 0.102 μm−1 (see Methods) and an increase in PSD at the lower frequency shoulder bordering the local minimum. The addition of a nucleus also impacts the slope of the high frequency region (βM) (Table 2 ) and skews it further to a value near 3 for all input β values.
Table 2. Impact of nuclear structure on βM values from representative SCO images
3.2 PSD-based characterization of simulated images containing epithelial structures
To examine more biologically relevant cell structures we threshold TPEF images of epithelial tissues, thereby isolating natural cell borders. We examine three layers of the epithelium, which have distinct morphologies: the differentiated (superficial), differentiating (para-basal), and undifferentiated (basal) layers (Fig. 3(a-c) ). We find that the average β exponents from fits to the PSD spectra of the original TPEF images of the superficial, para-basal and basal layers are 1.02 ± 0.18, 1.26 ± 0.20, and 2.01 ± 0.30, respectively. Differences between the exponents measured from the basal layers and the other layers are statistically different (p ≤ 0.003). To assess whether the differences in β exponents measured from these images originate from morphological variation between the cell layers and/or intracellular texture, TPEF binary masks generated by thresholding the images are analyzed. We find the average β exponents from fits to the superficial (2.30 ± 0.03), para-basal (2.48 ± 0.11), and basal (2.73 ± 0.13) layers of the masks are significantly different (p ≤ 0.032). This shift is related to cell morphology since there is an absence of intracellular texture in these images (Fig. 3(e-h)). For example, in our analysis of the basal layer images, the increased power of the PSD at frequencies corresponding to cell size, which is on the order of 0.1 μm−1, causes an increase in the slope of the power-law decay.
Fig. 3 (a-c) Representative original TPEF images from superficial, para-basal, and basal layers, respectively and (d) measured β exponents from corresponding PSD spectra and fits. Binary masks of cell borders (e-f) and corresponding PSD spectra and fits (h) demonstrating that a clear association of β exponents to morphology persists in the absence of high frequency spatial information. All images are 238 x 238 μm2.
To assess whether PSD-based techniques are sensitive to differences in subcellular texture in addition to cell morphology, fractals of varying β exponents are digitally applied to the intracellular regions of thresholded TPEF images from each tissue layer (Fig. 4 ). Figure 4 displays nomograms comparing the measured power-law exponents from fits to the PSDs of these simulated images to the input β exponents used in image generation. By digitally applying either a black background (high contrast, strong edge intensity) or the mean value of the image foreground (low contrast, weak edge intensity) to the background region, we are able to explore the effects of image feature edges on β exponent sensitivity. We find there is limited sensitivity of PSD-based methods to quantify intracellular texture in images containing a black background (Fig. 4(a-c), green lines) by observing that the measured power-law exponents of the fits varied marginally and non-linearly with input β exponents (Fig. 3). The average β-error, the average difference between βinput and βM for 0 < β < 2, for each tissue layer with a black background is 1.20 ± 0.74, 1.28 ± 0.75, and 1.63 ± 0.78 for the superficial, para-basal, and basal layers, respectively. By including the average value of the foreground in the background of the simulated images (Fig. 4(a-c); blue dashed lines), the sensitivity to input β values improves and average β-error is calculated to be 0.70 ± 0.32, 0.51 ± 0.32, and 0.69 ± 0.42 for the superficial, para-basal, and basal layers, respectively. However, in spite of this improvement, measured β values from the average value of the foreground in the background do not correspond linearly to input β values and are overestimated for β values lower than 3 and underestimated for β values greater than 3 for all tissue layers.
Fig. 4 Nomograms of measured β parameters for superficial (a), parabasal (b), and basal (c) epithelial layers of normal tissue. Plotted are input versus measured power-law exponents from linear fits of the radially-sampled PSDs for simulated images of varying input β parameters and background. Representative images with black background (green, solid lines), with the average value of the foreground in the background (blue, dashed lines) and images that have been clone stamped (red, dotted lines).
In order to obtain reliable sensitivity to a range of β parameters, a method to improve PSD sensitivity to input β values in the presence of cell features is necessary. Sensitivity to the correct β value depends on input β values, cell size, presence of a nucleus, and image background, which suggests that image features other than intracellular texture are contributing to the measurement of accurate power-law decay values. To recover β values that represent only the fractal nature of the intracellular signal, a digital object cloning (DOC) method is applied. DOC improves the accuracy of PSD-based characterization of intracellular structure. With DOC, fractal variation from simulated images can be recovered for β exponents that are less than the measured β exponents from the respective TPEF image masks (Fig. 4(a-c); red dotted lines). This limitation is tissue layer specific according to the measured β values from the image masks (Fig. 3). With DOC, relative average β-error decreases significantly to 0.06 ± 0.04, 0.05 ± 0.03, and 0.04 ± 0.04 for the superficial, para-basal, and basal layers, respectively.
3.3 DOC improves the accuracy of PSD-based characterization of intracellular texture allowing separation of normal and pre-cancerous tissues based on subcellular structure
After isolating cell features via thresholding, the DOC technique is applied to the original TPEF images from normal epithelial tissues. After DOC, the superficial, para-basal and basal layers (Fig. 5(a) ) have respective β values of 0.67 ± 0.13, 0. 69 ± 0.19, and 0.82 ± 0.14. β values between the layers are not statistically different (p ≥ 0.995).
Fig. 5 (a) Representative original and DOC-corrected TPEF images of superficial, para-basal, and basal layers of epithelial tissues made with healthy human foreskin keratinocytes (HFK) and HPV-transfected keratinocytes with the average β values for 5 different fields for each group displayed under each image with standard deviations. (b) β values from each HFK (green) and HPV (red) tissue sorted by tissue layer, before and after DOC.
DOC and subsequent PSD-analysis is applied to HPV-immortalized epithelial tissues (Fig. 5(a)). β values of original HPV tissue layers, shown in Fig. 5(a), are respectively 2.04 ± 0.20, 2.05 ± 0.08, and 1.98 ± 0.20 for the superficial, para-basal, and basal layers. After DOC, these values are 1.72 ± 0.19, 1.47 ± 0.15, and 1.53 ± 0.20 and are not statistically different from each other (p ≥ 0.227).
To visualize the effects of DOC across the tissue groups and layers, we plot the measured β values from TPEF images of each layer of five normal (Fig. 5(b), green) and five HPV (Fig. 5(b), red) tissues before and after DOC is applied. Before DOC is applied, there is not a significant difference between β values measured from HPV tissues layers and β values measured from the basal layer of normal tissue (p ≥ 0.963), which is most likely due to similarities in morphology and cell size. Remarkably, DOC corrects for tissue layer-dependent changes in β values for normal tissues by lessening the scale-dependent artifacts due to feature edges as well as image background. In this way, DOC enables targeted quantification of the fractal nature of subcellular structure, which is significantly different between the normal and the HPV-infected tissues layers after DOC is applied (p < 0.001). Although we find significant difference between tissue groups with a relatively small sample size, additional studies with larger samples sizes may provide insight into the diagnostic significance. In this way, DOC in combination with PSD-based methods shows great utility for biomedical image analysis, specifically with the goal of correlating subcellular structure to tissue health and disease status.
4. Discussion
This study demonstrates that Fourier-based analysis, relying on the PSD, is well suited to the characterization of high-resolution biomedical images for the assessment of intracellular structure. We specifically show that quantification of the β exponent of fractal biological features, such as mitochondrial organization, can be employed as an indicator of health and disease status. Assessment of mitochondrial structure may help elucidate mechanisms of normal and diseased tissue development, specifically in the context of cell metabolism. However, biomedical images generally contain a mix of fractal organization (e.g., intracellular morphology) and scale-dependent features (e.g., quasi-circular cellular and nuclear borders). The presence of nuclear and cell boundaries introduces edge-effect errors in PSD-based estimates of fractal character, which is clear when comparing pure fractal images (Fig. 1) with fractal images confined to circular shapes (Fig. 2). The images containing both scale-dependent and scale-invariant features require further processing to enable accurate characterization of the intracellular fractal patterns (Fig. 4). After standard thresholding-based processing of these images, fractal patterns with power-law decays greater than the power-law decay of the fractal pattern’s confining borders cannot be accurately measured due to edge artifacts. However, an additional DOC step allows reasonably accurate recovery of subcellular fractal organizational features.
4.1 DOC enables direct comparison of intracellular scale invariance among images with different large-scale features
Boundary effects due to overall cellular and nuclear shape, must be minimized in order to accurately quantify the fractal nature of intracellular structures. Decreasing the change in intensity at the edges of the cells by thresholding and applying the average intensity value of the foreground to the background (2-fold decrease in β-error) or thresholding and DOC (26-fold decrease in β-error) allows for more accurate PSD-based fractal characterization for β < 2. Despite the application of DOC pre-processing techniques to overcome edge artifacts, there is a drop-off in sensitivity to scale-invariant patterns characterized by power-law exponents greater than 2 compared to pure scale-invariant patterns (Fig. 4). Specifically, this occurs because our DOC method becomes less efficient at weakening edges within progressively more clustered images (β > 2), which have large deviations in local average intensity and are thus, more susceptible to DOC-induced edge artifacts (Fig. 1(e and f)).
Fortunately, the variation of intracellular organization observed from our analysis of experimentally-acquired TPEF images of engineered epithelial tissue occurs over a range of β values where PSD methods are highly accurate following DOC (β ≤ 2). In this regime, DOC reduces the contribution from the non-fractal features, allowing direct comparison of subcellular scale-invariance between images that contain large-scale, non-fractal variation. This is useful when evaluating different epithelial tissue layers that have varying characteristic cell sizes.
4.2 Comparison with other techniques
It is particularly advantageous to use DOC pre-processing to reduce edge effects because DOC does not sacrifice signal resolution, as with, for example, Fourier ‘windowing’ techniques [35]. Another Fourier-based technique developed with the goal of circumventing the limitation of Fourier analysis to rectangular images has been applied to the fractal analysis of nuclear chromatin distribution in benign and malignant breast cells [9]. This method fills the non-fractal region surrounding nuclear features using an iterative algorithm that terminates when the background region has similar statistical properties to the inner fractal region. For this method, it is difficult to obtain the background properties similar to the foreground, with 10% of samples unable to converge to the appropriate solution. Furthermore, the accuracy of the technique was evaluated using a single β value (β = 1.5), and it is unknown whether the method can perform as well with other β values.
Another method that is used frequently to determine fractal character is the box-counting approach [36, 37]. The box counting technique assesses the area that a binary signal covers at different scales. The method is performed by covering a thresholded image with a mesh of identical squares and counting the number of squares, N, that contain part of the image. This count is repeated for increasingly small squares (of size L) within the mesh. In this approach, a log plot N versus L is employed to investigate the scaling behavior of the putative fractal. This method is popular for its computational simplicity and for its ability to assess irregularly shaped images. However grid effects are known to interfere with computing the fractal dimension [36]. Further, box-counting methods often require total signal binarization using a simple threshold that must capture variations in fractal texture. It can be challenging to determine a single threshold value across a large data set [37]. Finally, fractal images with β in the range between 2 - 4 correspond to fractional Brownian functions for which box-counting methods are known to be inapplicable [28].
4.3 Diagnostic utility of combined DOC and PSD methods
Given the advantages of DOC for PSD analysis of intracellular texture in our simulations, we apply DOC to our original TPEF images. Using DOC, we find that PSD assessment of TPEF images has great promise for identifying diseased and healthy cells based on the ability to characterize the fractal nature of intracellular structures (Fig. 5(b)). When we compare the layers within normal tissues without DOC, we find significant β value differences between the healthy tissue layers. The β values measured from the HPV tissues do not change with tissue depth, but are statistically similar to the β values measured from the healthy basal tissue layer (p ≥ 0.963). Our simulation work suggests that this difference primarily arises from alterations in cell morphology, which is further supported by the similarities in morphology and β values between the basal layer of the normal tissues and all HPV tissue layers (Fig. 3). By assessing TPEF data after DOC, we find that that tissue disease status is related to the corresponding power-law exponent, with all normal tissue layers having significantly lower exponent values than HPV tissues (p < 0.001). Thus, the β value has potential to be used as a unique and valuable indicator of disease.
In the absence of edge effects, the β value represents the character of the fractal component of the image, which, in the context of this analysis, physically represents the organization of the mitochondria. We target mitochondrial organization by assessing images of NADH autofluorescence and by our choice of a fitting range that represents image structures on the size order of mitochondrial networks (< 10μm) [14]. We find an increase in β values measured from the analysis of precancerous (HPV) tissues, which reflects more correlated mitochondrial networks, as compared to the mitochondrial networks from normal tissues [38].
4.4 Disease progression and cell metabolism are related to mitochondrial autofluorescence patterns
Quantitatively assessing the nature of cellular mitochondrial networks and interpreting the physical meaning of the results has the potential to provide unique insight into the cell’s ability to satisfy its metabolic demands. Precancerous tissues are known to prefer glycolysis to oxidative phosphorylation, which we hypothesize is related to the increase in measured β values from our study [39]. Others have revealed that mitochondrial structure and internal organization is altered in response to the metabolic switch from glycolysis to oxidative phosphorylation [25]. From analysis of mitochondrial morphology using electron microscopy and fluorescence dyes, it was shown that a switch to oxidative phosphorylation is accompanied by higher order networking and an extension of the mitochondria into a more reticular form throughout the cell [25]. This is thought to be important in generating ATP at all parts of the cells. During glycolysis, the mitochondria were observed to cluster around the peri-nuclear area and were generally thicker. Thus, the increase of the β values throughout the layers of the pre-cancerous tissues could potentially provide insight on the altered metabolic demands of these cells towards a more glucose-oriented metabolism. This quantitative information provides insight into cellular metabolism that otherwise could not be attained with conventional analytical techniques and may have future diagnostic implications for automated real-time tissue characterization.
In addition to diagnostic implications, PSD-based quantification of biomedical images can help characterize complex biological patterns that may be essential to our understanding of the mechanisms behind differentiation, metabolism, and disease progression. Understanding the sensitivities and susceptibilities of these techniques to certain β value ranges and to image edges and feature morphology is essential in order to accurately assess images. Furthermore, in combination with non-invasive imaging techniques, this analysis enables quantification of subcellular patterns in live tissues without the use of added chemicals or contrast agents. The results of our analysis of autofluorescence images demonstrates that disease progression is related to the patterns intrinsically emanating from mitochondria, and that quantitative mitochondrial analysis may be useful for diagnostic purposes.
5. Conclusion
Fractals have become valuable descriptors of natural processes. Thus, there is a need for accessible analytical techniques whose limitations are fully understood, and at the same time, have the ability to accurately quantify scale-invariance found in nature. We find that the accuracy of the PSD-based techniques to characterize changes in fractal features within epithelial tissues is limited by the presence of nuclear and cell borders and image background signal. We develop a DOC technique that corrects for measurement inaccuracies induced by these edge effects, which results in a 26-fold decrease in measurement error and enables us to more accurately characterize naturally occurring mitochondrial fluorescence. By assessing the organization of mitochondrial networks, we are able to differentiate normal and precancerous epithelial tissues. We infer that these differences in organization are due to alterations in cell metabolic preferences. This technique could potentially be employed to correct other medical images with similar scale-dependent limitations, such as tissue histology or x-rays. Together, DOC and the PSD method have the ability to quantify the organization of a wide range of irregularly-shaped biomedical images and, more importantly, to isolate features of possible diagnostic significance. This could have applications in the development of automated disease detection software and computer vision.
Acknowledgments
This work was supported by American Cancer Society Research Scholar Grant RSG-09-174-01-CCE, NIBIB/NIH Grant R01 EB007542, and NIAMS/NIH Grant F32 AR061933. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
References and links
1. K. Doi, “Computer-aided diagnosis in medical imaging: historical review, current status and future potential,” Comput. Med. Imaging Graph. 31(4-5), 198–211 (2007). [CrossRef] [PubMed]
2. G. Dougherty and G. M. Henebry, “Fractal signature and lacunarity in the measurement of the texture of trabecular bone in clinical CT images,” Med. Eng. Phys. 23(6), 369–380 (2001). [CrossRef] [PubMed]
3. H. Gothwal, S. Kedawat, and R. Kumar, “Cardiac arrhythmias detection in an ECG beat signal using fast Fourier transform and artificial neural network,” J. Biomed. Sci. Eng. 4(04), 289–296 (2011). [CrossRef]
4. D. L. Turcotte, Fractals and Chaos in Geology and Geophysics (Cambridge Univ. Press, 1997).
5. P. Meakin, Fractals, Scaling, and Growth Far from Equilibrium (Cambridge University Press, 1998).
6. B. Mandelbrot, The Fractal Geometry of Nature (W.H. Freeman and Company, 2000).
7. H. S. Wu, “Fractal strain distribution and its implications for cross section balancing,” J. Struct. Geol. 15(12), 1497–1507 (1993). [CrossRef]
8. Y. Gazit, D. A. Berk, M. Leunig, L. T. Baxter, and R. K. Jain, “Scale-invariant behavior and vascular network formation in normal and tumor tissue,” Phys. Rev. Lett. 75(12), 2428–2431 (1995). [CrossRef] [PubMed]
9. A. J. Einstein, H. S. Wu, and J. Gil, “Self-Affinity and lacunarity of chromatin texture in benign and malignant breast epithelial cell nuclei,” Phys. Rev. Lett. 80(2), 397–400 (1998). [CrossRef]
10. S. A. Kartazayeva, X. Ni, and R. R. Alfano, “Backscattering target detection in a turbid medium by use of circularly and linearly polarized light,” Opt. Lett. 30(10), 1168–1170 (2005). [CrossRef] [PubMed]
11. J. M. Schmitt and G. Kumar, “Turbulent nature of refractive-index variations in biological tissue,” Opt. Lett. 21(16), 1310–1312 (1996). [CrossRef] [PubMed]
12. A. C. Sullivan, J. P. Hunt, and A. L. Oldenburg, “Fractal analysis for classification of breast carcinoma in optical coherence tomography,” J. Biomed. Opt. 16(6), 066010 (2011). [CrossRef] [PubMed]
13. M. Hunter, V. Backman, G. Popescu, M. Kalashnikov, C. W. Boone, A. Wax, V. Gopal, K. Badizadegan, G. D. Stoner, and M. S. Feld, “Tissue self-affinity and polarized light scattering in the born approximation: a new model for precancer detection,” Phys. Rev. Lett. 97(13), 138102 (2006). [CrossRef] [PubMed]
14. J. M. Levitt, M. Hunter, C. Mujat, M. McLaughlin-Drubin, K. Münger, and I. Georgakoudi, “Diagnostic cellular organization features extracted from autofluorescence images,” Opt. Lett. 32(22), 3305–3307 (2007). [CrossRef] [PubMed]
15. J. D. Rogers, I. R. Capoğlu, and V. Backman, “Nonscalar elastic light scattering from continuous random media in the Born approximation,” Opt. Lett. 34(12), 1891–1893 (2009). [CrossRef] [PubMed]
16. K. J. Chalut, J. H. Ostrander, M. G. Giacomelli, and A. Wax, “Light scattering measurements of subcellular structure provide noninvasive early detection of chemotherapy-induced apoptosis,” Cancer Res. 69(3), 1199–1204 (2009). [CrossRef] [PubMed]
17. M. Moscoso, J. B. Keller, and G. Papanicolaou, “Depolarization and blurring of optical images by biological tissue,” J. Opt. Soc. Am. A 18(4), 948–960 (2001). [CrossRef] [PubMed]
18. M. Bartek, X. Wang, W. Wells, K. D. Paulsen, and B. W. Pogue, “Estimation of subcellular particle size histograms with electron microscopy for prediction of optical scattering in breast tissue,” J. Biomed. Opt. 11(6), 064007 (2006). [CrossRef] [PubMed]
19. B. Chance, P. Cohen, F. Jobsis, and B. Schoener, “Intracellular oxidation-reduction states in vivo,” Science 137(3529), 499–508 (1962). [CrossRef] [PubMed]
20. B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts, and P. Walter, “Energy conversion: Mitochondria and Chloroplasts,” in Molecular Biology of the Cell (Garland Science, 2002) http://www.ncbi.nlm.nih.gov/books/NBK21063/.
21. C. R. Hackenbrock, “Ultrastructural bases for metabolically linked mechanical activity in mitochondria. II. Electron transport-linked ultrastructural transformations in mitochondria,” J. Cell Biol. 37(2), 345–369 (1968). [CrossRef] [PubMed]
22. C. R. Hackenbrock, T. G. Rehn, E. C. Weinbach, and J. J. Lemasters, “Oxidative phosphorylation and ultrastructural transformation in mitochondria in the intact ascites tumor cell,” J. Cell Biol. 30, 269–297 (1966). [CrossRef] [PubMed]
23. C. R. Hackenbrock, “Ultrastructural bases for metabolically linked mechanical activity in mitochondria. I. Reversible ultrastructural changes with change in metabolic steady state in isolated liver mitochondria,” J. Cell Biol. 51, 123–137 (1971).
24. H. Mortiboys, K. J. Thomas, W. J. Koopman, S. Klaffke, P. Abou-Sleiman, S. Olpin, N. W. Wood, P. H. Willems, J. A. Smeitink, M. R. Cookson, and O. Bandmann, “Mitochondrial function and morphology are impaired in parkin-mutant fibroblasts,” Ann. Neurol. 64(5), 555–565 (2008). [CrossRef] [PubMed]
25. R. Rossignol, R. Gilkerson, R. Aggeler, K. Yamagata, S. J. Remington, and R. A. Capaldi, “Energy substrate modulates mitochondrial structure and oxidative capacity in cancer cells,” Cancer Res. 64(3), 985–993 (2004). [CrossRef] [PubMed]
26. J. M. Levitt, M. E. McLaughlin-Drubin, K. Münger, and I. Georgakoudi, “Automated biochemical, morphological, and organizational assessment of precancerous changes from endogenous two-photon fluorescence images,” PLoS ONE 6(9), e24765 (2011). [CrossRef] [PubMed]
27. T. H. Wilson, “Fractal strain distribution and its implications for cross section balancing: further discussion,” J. Struct. Geol. 19(1), 129–132 (1997). [CrossRef]
28. R. F. Voss, in Fundamental Algorithms for Computer Graphics, edited by R. A. Earnshaw (Springer-Verlag, Berlin, 1985).
29. C. Meyers, T. J. Mayer, and M. A. Ozbun, “Synthesis of infectious human papillomavirus type 18 in differentiating epithelium transfected with viral DNA,” J. Virol. 71(10), 7381–7386 (1997). [PubMed]
30. W. R. Zipfel, R. M. Williams, R. Christie, A. Y. Nikitin, B. T. Hyman, and W. W. Webb, “Live tissue intrinsic emission microscopy using multiphoton-excited native fluorescence and second harmonic generation,” Proc. Natl. Acad. Sci. U.S.A. 100(12), 7075–7080 (2003). [CrossRef] [PubMed]
31. W. L. Rice, D. L. Kaplan, and I. Georgakoudi, “Two-photon microscopy for non-invasive, quantitative monitoring of stem cell differentiation,” PLoS ONE 5(4), e10075 (2010). [CrossRef] [PubMed]
32. J. Xylas, A. Alt-Holland, J. Garlick, M. Hunter, and I. Georgakoudi, “Intrinsic optical biomarkers associated with the invasive potential of tumor cells in engineered tissue models,” Biomed. Opt. Express 1(5), 1387–1400 (2010). [CrossRef] [PubMed]
33. R. E. Blahut, Theory of Remote Image Formation (Cambridge University Press, 2004).
34. K. P. Quinn, E. Bellas, N. Fourligas, K. Lee, D. L. Kaplan, and I. Georgakoudi, “Characterization of metabolic changes associated with the functional development of 3D engineered tissues by non-invasive, dynamic measurement of individual cell redox ratios,” Biomaterials 33(21), 5341–5348 (2012). [CrossRef] [PubMed]
35. F. I. Harris, “On the use of windows for harmonic analysis with discrete Fourier transform,” Proc. IEEE 66(1), 51–83 (1978). [CrossRef]
36. R. Lopes and N. Betrouni, “Fractal and multifractal analysis: A review,” Med. Image Anal. 13(4), 634–649 (2009). [CrossRef] [PubMed]
37. F. Normant and C. Tricot, “Method for evaluating the fractal dimension of curves using convex hulls,” Phys. Rev. A 43(12), 6518–6525 (1991). [CrossRef] [PubMed]
38. D. L. Turcotte, “Fractals in petrology,” Lithos 65(3-4), 261–271 (2002). [CrossRef]
39. O. Warburg, F. Wind, and E. Negelein, “The metabolism of tumor in the body,” J. Gen. Physiol. 8(6), 519–530 (1927). [CrossRef] [PubMed]
