1.

Introduction

The challenge of characterizing the seemingly random geometric relationships between the visible biological features of tissues can be simplified by examining the Fourier spatial frequency (FSF) spectrum of the image.¹ Drug delivery to tumors is well known to be chaotic and limited, partly from the perfusion limitations of dysfunctional neovasculature, but also because of the microscopic regional variations in composition which occur.^2,3 However, while the phenomenological features of drug delivery limitations are known, understanding ways to mitigate the problem remains an ongoing research challenge. The characterization of this seemingly random microheterogeneity is an issue that would benefit significantly from automated methods to quantify delivery efficacy. With this as a goal, in this study, the features of tumor sections were analyzed for anisotropic structures and how this relates to flow kinetics, by examining the FSF spectrum of the image features. This letter presents the conceptual development of how structure could be used to automatically classify transport in systems of well-characterized tracers and tumor types.

The intensity map of any image can be decomposed into the superposition of a series of sinusoidal frequencies, with units of cycles per unit length, similarly to the way electrical signals can be decomposed into individual temporal frequencies. Recently, Russell et al.⁴ used FSF analysis to identify the characteristic spatial frequencies of dispersions of self-assembled biological macromolecules in different states of aggregation. Pu et al.¹ demonstrated the use of FSF to distinguish between histological images of normal tissue and early-stage cervical cancer. Reufer et al.⁵ in a related application developed a method to incorporate directionality into FSF analysis, for identification of anisotropy in nanoparticle dynamics. Heterogeneous distribution of exogenous agents to solid tumor is thought to be associated with local tumor stromal heterogeneity: the relationship between nonuniform porosity from stromal growth and channel formation and nonuniform distribution of therapeutic or diagnostic agents has been observed both by histology,⁶ and in vivo MRI.⁷ The ultrastructure of the tumor stroma has been identified as a regulator of this phenomenon and is well known to be a major barrier in certain tumors, such as pancreatic cancer.⁸ Regions of invasive cell growth, disordered collagen, or high glocosaminoglycan content also can create a restricted geometry that defines a tortuous flow environment.⁹

In this study, we take advantage of the fact that in fluorescent tracer experiments, a region of poor tracer distribution demonstrates evidence of a barrier to transport. Using intensity maps of fluorescent tracers, we identify regions of marked fluorescence intensity anisotropy and then apply a directional averaging of the FSF over the intensity map of a conventionally stained histological section of that region. When a standard hematoxylin and eosin (H&E) stain is used, the eosin component, which stains protein, forms a pattern in the interstices between cells, i.e., an image of the protein distribution in the stroma. The characteristic frequency spectrum recovered from the intensity of the eosin staining is a quantitative record of how finely or coarsely spaced the structural features revealed in the image are. In this letter, we show that the characteristic spatial frequency spectrum of intratumor protein associated with directional anisotropy of tracer distribution is a unique identifier of transport barriers and establishs a systematic method of quantifying the relationship between tumor structure and drug distribution.

2.

Methods and Materials

Images were from cut sections of xenograft tumors grown in nude mice. All animal procedures were approved by the Dartmouth Institutional Animal Care and Use Committee (IACUC). A431 cells (ATCC, Manassas, Virginia) were implanted subcutaneously into the right flank via an injection of $1\times {10}^6$ cells in 50% Matrigel (BD Biosciences, San Jose, California) and 50% Dubelco’s modified essential media (DMEM) [10% fetal bovine serum (FBS), Pen/Strep] media, and grown for 8 to 14 days until reaching a volume of $\sim 200{\mathrm{mm}}^3$, then an injection of a $1:1$ mixture of 0.1 nmol of a epidermal growth factor receptor (anti-EGFR) Affibody (Affibody, Solna, Sweden) labeled with IRDye 800CW and 0.1 nmol of a nontargeted imaging control Affibody labeled with IRDye 680RD were made. Animals were sacrificed 1-h postinjection, the tumors removed, imbedded in Tissue Tek® Optimum Cutting Temperature medium (Sakura Finetek USA Inc., Torrance, California) and flash frozen. Tissues were stored at $-80°\mathrm{C}$ until sectioning. The frozen tumors were cut into 10 μm sections using a Leica 1850 M cryostat and fluorescence images of the tissues were immediately obtained on the Odyssey Scanner at a resolution of 21 μm (LICOR Biosciences, Lincoln, Nebraska) for both targeted and untargeted tracers. After fluorescent imaging, the tissue sections were stained with H&E by the Research Pathology Services at Dartmouth. The H&E stained tissues were imaged at $100\times$.

A schematic representation of key elements of the method is shown in Fig. 1. First, the intensity map of the nontargeted fluorescent signal was subtracted from that of the targeted Affibody, to normalize for regions of nonspecific tracer pooling. Then, the normalized fluorescent images were segmented into 256 intensity levels [Fig. 1(a)] and then aligned and overlaid upon the H&E images prepared from the same section [Figs. 1(b) and 1(c)]. Regions of interest were identified by visual inspection of the fluorescent contour maps [Fig. 1(d)], selecting regions showing a close spacing of contours. Axes parallel and perpendicular to the concentration gradient were identified [Figs. 1(e) and 1(f)] to apply directional Fourier analysis.

Fig. 1

Schematic illustration of the identification of contours and directional axes. (a) Scans of fluorescence from affibody-IRDye800CW were taken from unwashed frozen tumor sections and then segmented into 256 intensity levels for contouring. (b) Sections are then washed and stained with hematoxylin and eosin (H&E). (c) The segmented fluorescent images (green) were aligned with and overlaid upon the H&E sections (purple), to form full-section contour maps. (d) Detailed contour map of white box in C. (e) Focal regions of anisotropy were identified by steep contour gradients. The primary direction of analysis was perpendicular to the gradient (solid arrows) and was contrasted with the direction of steepest descent (dotted arrows). Arrows indicate the direction of analysis only, not direction of tracer flow. (f) Axes were overlaid upon the corresponding H&E image for FSF analysis. Note that the tumor in these images shows evidence of tissue loss from the presence of a necrotic core, a common feature of this xenograft. It is significant that since the H&E stain is prepared subsequently to the fluorescent image, from the same tissue, the interior limit of the fluorescent signal is not a result of tissue loss during the staining process but represents limited tumor penetration.

These regions of the H&E images were evaluated for directional anisotropy by azimuthally averaged frequency power spectrum analysis, using a variation of the method introduced in Ref. 5. The method of generating directional FSF spectra is as follows: after identification of the region of anisotropic tracer concentration, a square patch of the tissue section, centered on the point at which the two axes of anisotropy intersect [Fig. 2(a)], is selected and rotated so that the axis of interest was at $-45\mathrm{deg}$ [Fig. 2(b)]. The three-color image was converted to an intensity map, as shown in Fig. 2(c), and the intensity of the cell nuclei was set to zero. This had a dual effect: first, it reduced high-frequency contributions from nuclear structures that were unlikely to have an effect on tracer transport; and second, it ensured that spatial frequencies corresponding to gross spacing between cell nuclei would be below 5 cycles/cm, allowing elimination of these with a low-pass filter. The two-dimensional (2-D) FT of the rotated intensity map was then generated by the built-in MATLAB function FFT2.

Fig. 2

Schematic illustration of the method of generating directional Fourier spectra. (a) A square patch of the H&E image was selected centered at the point of intersection of axes of anisotropy. (b) The patch was rotated so that the principal axis of analysis is at $-45\mathrm{deg}$. (c) The H&E image was converted to an intensity map, and cell nuclei were assigned an intensity of zero. (d) The 2-D FT of the rotated intensity map is generated. The resulting function was numerically integrated at steps of 1 deg over an arc of 7 deg centered on the principal axis of analysis, as shown by the short solid arrow, and shaded region R. Spatial frequencies characteristic to the image in the region of analysis were identified by numerically integrating the remaining 352-deg arc of the full circle, again at 1 deg steps, as shown by the long solid arrow, and shaded region W. The frequency power spectrum of the secondary axis of analysis was generated by a similar numerical integration, as shown by the short dashed arrow and shaded region, and the image-specific spatial frequencies found by numerical integration of the remainder, as shown by the long dashed arrow.

3.

Theory

The power spectrum of an intensity map $M(u,v)$, is written as ${|F(u,v)|}^2=R^2(u,v)+I^2(u,v)$,¹⁰ where $R(u,v)$ and $I(u,v)$ are the real and imaginary parts of the FT, respectively. Following the method established in Ref. 4, we decompose the power spectrum into three components, so that ${|F(u,v)|}^2=T(u,v)+S(u,v)+B(u,v)$, where $T(u,v)$ is the characteristic frequency that defines similarity between image regions, $S(u,v)$ is the contribution of the spectrum that is specific to the region, and $B(u,v)$ is a noise component that comes from instrumentation artifacts and the analytical limitations of the finite Fourier transform. To extract $T(u,v)$ in the direction of interest, the azimuthally averaged power spectrum of a 7-deg arc of a sub-region $R$ of the 2-D FT centered around the principal axis of analysis was generated by numerically integrating the signal at 1 deg steps. The background signal was generated by integrating over the remaining 352 deg of the 2-D FT, sub-region $W$, in Fig. 2(d), again by numerical integration at steps of 1 deg. Because of the high variability in FT traces over the large range of $W$, the expected contribution of $T_W$ and $S_W$ is zero, and by definition $B_R=B_W$, giving

4.

Results

In Figs. 3(a) and 3(b), individual FSF spectra from the principal and secondary axes of analysis of a representative image are plotted on the same axis, and the geometric mean of the sample is overlaid on them. The natural variation expected from a biological sample can be seen, as well as a striking correlation between the individual traces and the mean. In Fig. 3(c), the power spectrum of the two directional spectra is plotted on the same axis, showing spatial frequencies unique to each direction, with very little overlap. The method was applied to several regions of interest derived from contour maps, without regard to any apparent similarity of the underlying H&E staining. Figure 4 shows representative images with similar structures of tracer distribution anisotropy, which overlie regions of tumor stroma that have a markedly different visual appearance. The specific frequencies that identify stromal structural elements are not identical between different regions, but the characteristic directional spectra in the direction of flow can be seen to be distinct from those in the direction perpendicular to flow in each case.

Fig. 3

Representative directional frequency spectra generated from intensity maps as described in the text and shown in Fig. 2. (a) The spectra in the direction parallel to the fluorescent tracer flow taken from eight points along the path shown in Fig. 1(f) are plotted on the same axis in gray. The geometric mean of the traces, defined as $\prod F_i$ ($1/i$), where $i=8$ in this case), is shown by the solid line with open circles. (b) The same presentation as A, but in the direction perpendicular to the local fluorescence gradient/flow, and the geometric mean shown by the solid line with solid circles. The geometric mean is plotted against a second $y$-axis. (c) The averaged power spectrum, being the inner product of all traces, and in the same markers as the geometric mean in (a) and (b). The directional FSF in the direction of steepest descent shows very poor correlation at frequencies below $2\text{cycles}/\mu \text{m}$, whereas the directional FSF parallel to the intensity contour shows low amplitude at frequencies above $2\text{cycles}/\mu \text{m}$.

Fig. 4

Sample contour map images, and directional fourier spatial frequency analysis. (a) Contour map of fluorescent tracer distribution overlayed on an H&E image prepared from the same tissue section (A431). Principal direction of analysis is indicated by the white arrows (arrows indicate direction of analysis only, not direction of flow). (b) Fourier spatial frequency spectra in the principal direction of interest (open triangles) and the perpendicular direction (solid circles). (c) A cumulative sum of the frequencies in (b). The separation of the curves shows the predominance of high frequencies in the direction of interest, relative to the perpendicular. (d) Another contour map of the same tumor, with the principal direction shown by the white arrows. (e) Markers as in (b). The spectra are distinct from those in b but still show distinct characteristic spectra by direction. (f) Cumulative sum shows separation at low frequencies. Spatial frequencies have units of cycles per μm.

5.

Discussion

We show here that the direction of free tracer transport can be differentiated from that of impeded tracer transport by comparison of the different characteristic spatial frequencies of the underlying eosin stain pattern. The system we use, H&E histopathology, is one that is well established in evaluation of cancer in the clinic. The pink stain of eosin, typically associated with protein, provides the signal from which we extract the characteristic directional FSF spectrum. The evaluative potential of this method can be extended by developing a staining design strategy that minimizes extraneous information, for example, by eliminating the hematoxylin stain, or that adds information, such as by locating microvessels by staining for vascular components.

It should be noted that the individual frequencies contained in a power spectrum do not have a straightforward correlation to features in an image, except in the case of an image with a regularly repeating sinusoidal pattern. For instance, the presence of a cluster of many high spatial frequencies is not necessarily an indicator of a large number of small features. It would be a fundamental misapplication of the method to expect that the presence of specific frequencies is an indicator of high or low transport. The correlation between spatial frequencies and visual patterns forms a distinct field of optical research and is beyond the scope of this letter. An excellent introduction to the topic can be found in Ref. 11.

Although there is an intuitive sense of what “the same” means from a pathologist perspective, the question of how to define similarity and difference, is not specific to spectral analysis or even to the natural sciences and can be found underlying arguments as diverse as biological taxonomy and the literary theory of deconstruction.¹² For the purpose of FSF analysis, the key to identifying a characteristic spectrum in a particular direction is to differentiate the signal from a suitable background signal. The orthogonal spectrum, because it uses the same background, can always be distinguished from the principal direction if it has a distinct spectrum.

Because the fluorescent overlays show tracer distribution to be anisotropic, and there is some indication of this anisotropy in the H&E images even before analysis, the results presented here are a quantitative confirmation of what can be surmised. There are two principal experimental extensions of this method. First, as shown here, it provides a robust quantitative characterization of spatial anisotropy in systems where it can already be qualitatively detected. Directional FSF analysis can also be applied to imaging systems where a spatial anisotropy is expected but difficult to detect. Such would be the case in two-photon (2P) imaging of biological tissue, in which a second harmonic of the emission frequency is generated by collagen fibrils, and a variation in signal can be expected depending on the ordered state of the collagen [reviewed in Ref. (13)]. One distinct advantage of this application of 2P imaging is that it allows for real-time monitoring of drug distribution, while allowing for simultaneous capture of the signal used for FSF analysis. Likewise, it may be possible to monitor real-time destruction of transport barriers by activable nanoparticles¹⁴ designed for that purpose. A further extension of the method can be made to a 3-D case by stacking 2-D image slices.

In principal, the method of directional FSF analysis can be applied to any signal or imaging system. But there are two prerequisites for such a system. First, the signal to be analyzed must be available as a spatially distributed intensity map in which the intensity signal has some (known or unknown) relationship to the question at hand. Without this format, an informationally relevant FSF spectrum cannot be generated. Second, a suitable method must be identified for generating a background signal against which to isolate the signal of interest. This is especially important for biological applications, in which background noise can be on the same order of magnitude as the information.

In many tissues, structural anisotropy is strongly correlated with functional parameters. Diffusion tensor MRI can identify tracts of bundled neurons¹⁵ and muscle fibers,¹⁶ by the preferential diffusion of water along them, and automated methods have been explored.¹⁷ MRI tractography is possible because of the physics of the detection system, which allows for a large number of measurement directions (typically around 20) and the acquisition of a large number of data sets (typically thousands). An optical analogue of MRI tractography would be of significant benefit in the detection of anisotropic drug delivery in tumor. A key component will be an automated quantitative method of identifying similarities and differences between adjacent flow regions in optical maps, which we have demonstrated here. In future work, we hope to apply this method to images stained for microvasculature, to allow a discussion of fluid flow as it pertains to barriers to drug distribution. In addition, we will generate an automated scanning protocol that will examine steps along the track of the direction of interest to detect paths that have similar spectral characteristics that branch from or extend the initial path.