Significance: Tumor heterogeneity poses a challenge for the chemotherapeutic treatment of cancer. Tissue dynamics spectroscopy captures dynamic contrast and can capture the response of living tissue to applied therapeutics, but the current analysis averages over the complicated spatial response of living biopsy samples.
Aim: To develop tissue dynamics spectroscopic imaging (TDSI) to map the heterogeneous spatial response of tumor tissue to anticancer drugs.
Approach: TDSI is applied to tumor spheroids grown from cell lines and to ex vivo living esophageal biopsy samples. Doppler fluctuation spectroscopy is performed on a voxel basis to extract spatial maps of biodynamic biomarkers. Functional images and bivariate spatial maps are produced using a bivariate color merge to represent the spatial distribution of pairs of signed drug-response biodynamic biomarkers.
Results: We have mapped the spatial variability of drug responses within biopsies and have tracked sample-to-sample variability. Sample heterogeneity observed in the biodynamic maps is associated with histological heterogeneity observed using inverted selective-plane illumination microscopy.
Conclusion: We have demonstrated the utility of TDSI as a functional imaging method to measure tumor heterogeneity and its potential for use in drug-response profiling.
Tumor heterogeneity presents a challenge for the successful treatment of cancer using chemotherapeutics.1 For instance, genetic variability in tumors caused by clonal outgrowth of selected genotypes within a tumor may cause subsets of cells with genetic variations to be resistant even while the majority of the tumor responds to treatment. Selective pressure and genetic drift of the cancer cell population during treatment often lead to patient relapse and the emergence of broad chemoresistance and refractory disease.2–4 In addition to genetic heterogeneity, there is also spatial heterogeneity in tumor tissue arising from varying tissue constituents as well as varying microenvironments, including differences in extracellular matrix and connective tissues. The tumor microenvironment5,6 and epigenetic variations5,7–9 pose significant challenges to the selection of treatment based on genetic profiles. This has led, as an alternative, to phenotypic profiling10–12 of cancer tissue that captures the systemic response of cancer tissue to applied therapy. The challenge for phenotypic profiling of cancer tissue is the need to image intact microenvironments deep inside tissue, far from surface damage caused by surgical resection, and deep inside transport-limited regions that experience hypoxia, nutrient depletion, and metabolite build-up.
Optical coherence imaging (OCI)13,14 is a deep-tissue coherence-domain imaging approach based on digital holography15–17 that is a form of full-frame optical coherence tomography.18,19 Dynamic speckle in OCI images caused by dynamic light scattering from intracellular motions enables biodynamic imaging (BDI)20 to use intracellular dynamics as a unique form of image contrast. The changes in intracellular motions caused by applied therapeutics have been studied using tissue dynamics spectroscopy (TDS)21 to separate the effects of drugs across broad spectral bands and to capture specific signatures from different classes of drugs with different mechanisms of action.22 Preclinical trials of therapy responsivity assessment have been completed using TDS in spontaneous canine B-cell lymphoma and in ovarian xenografts.23,24 In a substantially different setting, assisted reproductive technology correlates the viability of cumulus-oocyte complexes with parameters from sample fluctuation power spectra.25
The methodology of TDS is usually applied to entire samples that can be as large as in volume (e.g., biopsies). However, intrasample variability in the TDS signatures poses a challenge for the prediction of patient response to therapy. While previous work using TDS has identified and characterized the different baseline conditions and drug responses in the “shell” and “core” areas of the samples,21,22,25 in that analysis, the boundary between the shell and core was arbitrarily defined. Some samples have a more complicated drug response structure than a simple “shell” and “core” model, as shown in Fig. 1, where the sample shows variation in both strength and pattern in its drug response in the two areas.
To address this problem, we introduce a functional imaging method called tissue dynamics spectroscopic imaging (TDSI) that evaluates sample response on a pixel level. In addition to a full-duration response map, the response can be segmented along the time axis to derive the time-lapse evolution of drug response, which can reveal the different rates at which a drug acts on each area. This methodology offers a quick, intuitive visualization of sample heterogeneity and drug effects. TDSI differs from dynamic-contrast OCT performed with off-axis holography,26 full-field OCT,27 or conventional OCT28–33 because its imaging contrast arises from shifts in intracellular dynamics caused by applied therapeutics rather than from baseline intracellular dynamics. In this way, TDSI is functional imaging that is specific to drug efficacy and displays the heterogeneous response of tissue to therapies.
Materials and Methods
Biological samples used in this paper include tumor spheroids grown from the DLD-1 intestinal adenocarcinoma cell line (ATTC, Manassas, Virginia) and human esophageal tumor biopsies. (IRB-approved IUCRO-0486 clinical trial. Written consent was obtained for each patient.) The multicellular DLD-1 spheroids mimic small avascular in vivo tumors by having an active growth zone on the outside of the spheroid and regions of apoptotic and necrotic cells progressively toward the core as nutrients and oxygen become rate-limiting for cell growth, but otherwise the spheroids are highly homogeneous. In contrast, the esophageal tumors are highly heterogeneous, with complex physiology and consist of different cell types. DLD-1 spheroids were grown in Corning U-bottom 96-well spheroid plates, and esophageal tissues were obtained by pinch biopsy from human patients. Tumor biopsies were collected and transported in chilled RPMI-1640 medium with HEPES buffer and cut into small pieces of 1 mm size or less. Both types of samples were immobilized in 1% low-gel-temperature agarose in the RPMI-1640 basal medium. Immobilized samples were overlaid with RPMI-1640 containing 10% heat-inactivated fetal calf serum (Atlanta Biologicals), penicillin (100 IU), and streptomycin ().
Sample heterogeneity observed using TDSI was verified with high-resolution three-dimensional (3-D) images obtained from inverted selective plane illumination microscopy (iSPIM). To prepare samples for iSPIM imaging, 10% neutral buffered formalin was injected into each well to fix the tissues. After being washed with PBS, the samples were stained with DRAQ5 (Biostatus, Ltd.) overnight with gentle agitation, then 80% ethanol-based Eosin Y (E4009, Sigma-Aldrich) for 4 h. Afterward, the agarose embedding the samples in the dish wells was removed, and samples were then washed with DI water three times, and PBS once followed by being immersed in X-CLARITY mounting solution (Logos Biosystems) for 15 min. Finally, the samples were fixed in the imaging chamber with silicone glue and immersed in X-CLARITY mounting solution for iSPIM imaging. After being imaged, the samples were processed for traditional H&E through the Tulane Medical School Histology Department. Four-micrometer-thick sections were cut and stained until each tissue was exhausted. The total processing time is about 16.5 h in total after fixation, including the PBS/di-water washing time. Imaging of each sample takes about 1 min average scanning time. Full-resolution image processing of each sample takes , including time for reconstruction, making pseudocolor images, and visualization.
Tissue Dynamics Spectroscopy
Speckle fluctuation dynamics of biological samples were measured and analyzed using a BDI system. The system optical configuration is a Mach–Zehnder interferometer, shown in Fig. 2. The light source is a Superlum S840-B-I-20 superluminescent diode with a center wavelength at 836.2 nm and a full power output of 22.9 mW with a bandwidth of 50 nm and a coherence length of . A QImaging EMC2 camera is used for image acquisition. Lenses L5 and L6 constitute a imaging system to an image plane (IP). The IP is subsequently transformed by lens L7 to the Fourier plane located at the CCD pixel array. A single sample image is obtained by performing a two-dimensional (2-D) FFT on the hologram captured by the CCD camera located on the Fourier domain, and one of the two conjugate images is stored as a array.21 Because of the use of long focal lengths (), the speckle size on the camera plane is , and the reconstructed object-plane point-spread function is spanned by three pixels.
Two data acquisition formats are used for data presented in this paper: a format containing 2048 frames captured at 25 fps, and a format with 500 frames at 25 fps and 50 frames at 0.5 fps. The second format has a smaller data storage requirement but requires a stitching algorithm to create a single spectrum.34 Within this paper, only DLD samples use the 2048 frame format, while the esophageal samples used the 500/50 frame format to assure comparability among spectra within a given study. A typical experiment has six time-frames of baseline measurements and 15 time-frames of drug response measurements looping repeatedly through 16 wells in a 96-well plate format. A single loop through all 16 samples takes about 40 min.
The fluctuation power spectrum from one pixel at position in the sample is calculated as the square of the FFT of the intensity time series:
The average spectrum of a region is calculated as34
Biodynamic Biomarkers and TDS Visualization
Biomarkers that evaluate sample preconditions, also called baseline conditions, before a treatment is applied include backscatter brightness, normalized standard deviation, spectrum knee frequency, and spectral slope (S).32 Biomarkers that evaluate drug responses include changes in these biomarkers after application of drugs as well as features extracted from spectrograms. This paper focuses mainly on the drug response and features from spectrograms.
Despite the differences in drug responses related to sample baseline conditions and drug mechanisms, a biodynamic drug spectrogram usually has one of a limited number of patterns. Spectroscopic masks (time-frequency filters) are designed to match the characteristics of the spectrograms, a few of which are shown in Fig. 3(a). The top three patterns G0, G1, and G2 form a set of “orthonormal” masks that are related to the broadband (in the sense of frequency components) pattern of a spectrogram, while the bottom three patterns FL, FM, and FH form another set of masks related to local response patterns. These frequency bands, and their associated frequency cutoffs, can be related to changes in intracellular motions and their related speeds.23,35 Although the individual biomarker values depend on the choice of cut-off frequencies, the orthonormal character of the masks provides a unique representation of the spectrogram, and this representation is used in pattern recognition algorithms. For each mask, a feature value is obtained by calculating the inner products of the spectrogram and the mask, i.e., projecting the spectrogram onto the mask,22 and the features of a spectrogram are represented by a vector of feature values.
After condensed-data-format files are generated (refer to Sec. S.1 in the Supplemental Material), a differential spectrogram for each TDSI pixel (referred to as a “microspectrogram”) is calculated asFig. 3(b). For the “G0” TDS image, the area in the blue circle has negative values, indicating broadband decrease in activity, which matches the differential spectrogram of the circled area of the same sample as shown in Fig 1(a). Similarly, the spectral response in the red circle area agrees with the positive values in the same area in the “G1” TDS image.
As discussed above, the sample in Fig. 3 has a large variation in drug response in terms of strength and spectral patterns, displayed by the distribution of values in the TDSI of Fig. 3(b). Both “G0” and “G1” images show a change in magnitude and sign from bottom left to top right (corresponding to changes in the strength of drug response, referred to as “intramask heterogeneity”), and the “G0” image has a different pattern than “G1,” where “G0” has strong negative values on the bottom right while “G1” has strong negative values in the upper left (corresponding to changes in pattern, referred to as “intermask heterogeneity”). To better visualize the variation, bivariate color images are introduced to produce a single visualization that captures drug responses across the entire sample, where each “variable” is a feature value of a mask. Feature values from two masks are a convenient way to illustrate drug-response heterogeneity within a sample, and the following discussion will focus on bivariate representations of drug response.
Bivariate color maps are used in cartography36,37 and medical imaging,38 and many studies have addressed how to choose proper colormaps for bivariate data visualization.39–41 We have elected to use the “Teuling3” colormap in the following visualizations, which is generated by linearly interpolating four colors at the four corners in the sRGB space plus a “whitening” core in the center.39,42 This color map has good color saturation, relatively equal visual impact, and a zero value appears as white, which is consistent with the diverging “blue-red” one-dimensional colormap used in our spectrograms and univariate TDS maps. Figure 3(c) shows a bivariate image of a human esophageal biopsy sample “merged” from the two univariate TDS maps in Fig. 3(b).
When a sample has areas with spectrograms that are the inverse of each other, they cancel each other out in an average over the full sample, resulting in a mild spectrogram and near-zero feature values. In this case, the weak average response belies the strong change in the intracellular dynamics and the fluctuation spectra, and this can lead to the misinterpretation that the sample does not respond to the drug. To address this problem, two new biomarkers that evaluate sample heterogeneity are added to the “traditional” average spectrogram-based biomarkers. The two biomarkers evaluate the “intramask” and “intermask” heterogeneity, respectively. To achieve a high signal-to-noise ratio, TDS images are first (re)generated with an averaging (instead of the standard ). The coarse eight-pixel averaging is used only to calculate the heterogeneity biomarker values of a sample at a scale of and above. Finer scale is not necessarily meaningful for quantifying spatial heterogeneity across a millimeter sample. The feature values are bounded to a range and then assigned “scores” ranging from 0 to 1 for both heterogeneity evaluations, calculated using the following steps:
1. Select a set of masks
2. For each mask , calculate and
3. For each mask pair , calculate and
4. The first biomarker denoting overall intramask heterogeneity is calculated as
5. And the second biomarker representing overall intermask heterogeneity is calculated as43 We use and , and the local masks FL, FM, and FH are used for heterogeneity scores.
Based on these definitions, the extreme values are achieved under these cases: when TDS images are completely uniform (totally homogeneous), and when TDS images have only two values of opposite signs in an equal number of pixels. When the TDS images are completely nonlinearly correlated to each other, then . In the opposite case when all values are completely correlated then . Examples will be provided in the next section to illustrate these two heterogeneity benchmarks.
Time-Lapse Drug Response Visualization
After a drug is added to a biopsy sample, the change in its intracellular dynamics is usually not immediate and depends on the drug mechanism of action, especially for drugs targeting DNA that produce slow cellular responses. Also, some parts of the sample may respond to a drug faster than the entire sample. TDSI allows us to study both the time delay and nonuniformity in drug action, which is called time-lapse TDSI. Instead of extracting feature values from full-length spectrograms, time-lapse TDSI uses responses within a small moving time “window” of the spectrogram. Examples are included in the following sections.
Inverted Selective Plane Illumination Microscopy
The iSPIM system is constructed around a commercially available diSPIM platform (ASI) and has been described in previous publications.44,45 In brief, two immersion objectives (CTO, ASI/Special Optics, ) are orthogonally mounted above the sample, with each at a 45-deg angle from the norm, enabling traditional sample preparation. Dual-view imaging is possible by alternating roles of the two objectives as illumination and detection, but only single-view was adopted for this paper. Volumetric images were obtained by moving the sample with respect to the objectives. To cover the whole area of the sample, multiple y strips were acquired with about 20% overlap between two adjacent strips. After imaging was completed, the images were first shifted and interpolated with custom MATLAB code to recover its 45-deg angle.46 Multiple paths were then stitched with Fiji plugin.47 The DRAQ5 and eosin images (D&E) were remapped to a composite RGB stack to simulate the traditional H&E colors.48 Finally, 3-D reconstruction of the D&E images was obtained from the alpha blending mode of 3-D viewer of Vaa3D.49,50 The mounting media is Xclarity with a refractive index of 1.45, that when combined with the NA of 0.4, gives a magnification of with a lateral resolution of and an axial resolution of .
A large number of esophageal biopsies display spatially heterogeneous responses to drugs. In Fig. 4(a), two biopsy samples that have large intramask heterogeneity are presented in univariate and bivariate forms. Sample “151208-6” has a weak response when averaged over the whole-sample spectrogram [Fig. 4(b) “global”], because the local areas “1” and “2” have strong but opposite responses that tend to cancel in the sample average.
An assortment of bivariate TDS images is shown in Fig. 5. Some samples have relatively uniform color in the “merged” map, indicating smaller variation in the biomarker values, while others have a rainbow-like smooth transition across the sample, which is related to high heterogeneity in the drug response.
There are roughly three types of heterogeneity, shown in Fig. 6 along with their and scores: (i) Type I are samples that have almost uniform responses under all masks. (ii) Type II samples have spatial variability but show similar patterns across different masks. (iii) Type III samples have TDS images with nonoverlapping strongly responding areas. Types I and III are the most common.
Time-lapse images offer an additional layer of understanding of the spatial evolution of drug effects or sample conditions. In Fig. 7 for sample “170317-9” treated with nocodazole, the blue response pattern (mid-frequency suppression and low- and high-frequency enhancement) grows stronger over time before saturation, which indicates that nocodazole’s suppression of microtubule polymerization begins at the outer periphery and slowly penetrates the core of the sample. As another example, the red area in the TDS image of sample “170606-15” becomes stronger until around 9 h, when the sample displays an overall suppression across the entire sample.
Comparison with Inverted Selective-Plane Illumination Microscopy
Given that BDI is a 3-D imaging technology that uses low-coherence light, the different drug response phenotypes revealed by TDSI may be related to different types of tissues in a certain region of a sample. A complementary 3-D imaging technique is iSPIM that produces microscopic images of 3-D slices with high lateral and axial resolutions, which allows us to distinguish features in the images. Therefore, by comparing TDSI with iSPIM, we can investigate whether the heterogeneity related to drug response variability from TDSI is also present in microscopic images, i.e., link dynamic information from functional imaging with the histology of biological tissues.
As an example, TDSI maps for sample 190801-15 are compared against its iSPIM images and H&E histology images in Fig. 8. This sample is a pinch biopsy from a patient who was resistant to neoadjuvant therapy. The biodynamic response is mapped for the tissue response to 5-fluorouracil. In the TDSI map, there is a distinct central region (purple) contrasted with the outer regions (orange). In the bivariate colormap, purple represents broadband excitation with a blue shift, indicating an activated response of the tissue to the 5-fu. The central region is likely to be naturally hypoxic, which can affect the mechanism of the drug. In a related study of esophageal cancer patients, a blueshift induced by 5-fu is representative of a beneficial response to the chemotherapy. In contrast, in the peripheral tissue, that is hyperoxic, the broadband excitation is accompanied by a redshift. Furthermore, the green region indicates overall suppressed activity (broad inhibition accompanied with a redshift of slower intracellular speeds). In comparison to the TDSI, in the iSPIM image the lower part matches the lower part of the histology image containing a large concentration of DNA. The upper part of the iSPIM image is cytoplasm or unstained tissues, matching the lack of nuclei in the upper part of the histology image, which potentially indicates collagenous connective tissues. Although the image orientations were not registered between the techniques, the images in Fig. 8 demonstrate spatial polarity in the tissue types. Future studies to compare TDSI with iSPIM would register sample orientation to identify TDS spectral signatures of connective tissue relative to epithelial tissue. The work presented here is proof-of-principle that spatial heterogeneity can be observed in both imaging modalities.
In this comparison, the advantage of iSPIM is the ability to acquire microscopic image with high resolution and detailed cellular-level structural information with acceptable speed. This is in contrast to TDSI which achieves only tissue-scale imaging resolution ( size). On the other hand, TDSI is a functional imaging method that is highly specific to the action of drugs on tissue dynamics, while iSPIM can only provide structural information without being specific to drug response. Furthermore, TDSI is nondestructive, allowing longitudinal studies that can span several days, including the ability to apply drugs and to clear them. Both techniques have the advantage of imaging into 3-D tissues, which is a key aspect of maintaining the tumor microenvironment during imaging.
BDI is a tool that is sensitive to intracellular dynamics and has been applied successfully to phenotypic profiling to predict patient outcomes. For instance, sample motility and dynamic biomarkers have been shown to be consistent and reliable indicators of pharmacodynamics effects. However, these biomarkers are usually calculated as whole-sample averages when used in classification and similarity analyses, overlooking intrasample heterogeneity. The introduction of TDSI solves this problem by evaluating the responses of subregions of the sample to reveal new information on the complex spatial structures in sample drug response, which is supported by evidence from other imaging techniques such as iSPIM and histology. Visualization of sample heterogeneity is facilitated with a bivariate color representation and is quantitatively characterized by and scores. This imaging method is further extended to generate whole-sample time-lapse TDSI maps, providing a method to monitor drug mechanisms.
In addition to visualizing sample heterogeneity, TDSI maps can provide additional information and improve classification accuracy when evaluating anti-cancer drug effectiveness on a patient level. For samples with regions of opposite responses, the whole-sample average spectrogram may suggest a mild response to the drug, making the sample and the patient appear to be less sensitive to the treatment. The proposed solution here is to introduce additional biomarkers that characterize regional drug responses. As an example, the sample shown in Fig. 3(b) can be split into two regions based on the sign of G1 biomarker values (which can be related to the sample’s heterogeneous structure), and the feature values of these two regions can be calculated. The set of feature values that capture both sample average response as well as regional response would provide a more comprehensive assessment of the patient.
TDSI can be extended for further imaging and analysis applications. Since BDI is a 3-D imaging technique and achieves depth selection with coherence gating, a volumetric TDS image can be generated by scanning different slices of a sample. Also, time-lapse TDS analysis is a quantitative approach to visualize drug-action time dependence. Features such as delay and distribution related to pharmacokinetics and pharmacodynamics (PK/PD) could be obtained from time-lapse images to provide insight into processes such as dose-response relationships.
Challenges to TDSI include sample immobilization and multiple light scattering. TDSI evaluates the drug spectral response on the pixel level and requires that the same part of the sample is imaged throughout the experiment. This requires the sample to maintain the same lateral and axial positions. In addition, multiple light scattering induces aberrations of the image and reduces the signal-to-noise ratio, which makes TDSI more effective at shallower depths.
TDSI is an important extension to the current suite of BDI modalities (OCI, MCI, and TDS). MCI maps are simple and intuitive functional images that visualize sample motility and have revealed the contrast between a viable shell and a necrotic core for rat tumor spheroids.26 TDSI, by comparison, generates a set of more detailed functional maps that complement MCI because the critical frequencies in spectral masks used in TDSI are related to specific types of intracellular components and motions, offering a comprehensive view of changes occurring in the sample. TDSI is a versatile functional imaging method that could provide new information for drug-response profiling and has the potential for improving predictions of response to therapy and drug screening.
David Nolte and John Turek have a financial interest in Animated Dynamics, Inc., a company that is seeking to commercialize biodynamic imaging for chemotherapy selection.
This work was supported by the U.S. National Science Foundation (NSF) and National Institutes of Health (NIH), Grant Nos. NSF 1911357-CBET and NIH NIBIB 1RO1EB016582, and by the Purdue University Cancer Center.
Zhe Li received his PhD in physics from Purdue University in 2020. He studied intracellular dynamics of living tissues using Doppler fluctuation spectroscopy under the direction of Dr. David Nolte. He currently works on model verification and validation products at MathWorks, Natick, Massachusetts, USA.
Bihe Hu received her BS and MS degrees in biomedical engineering from Huazhong University of Science and Technology, where her research interest was in the area of 3-D whole-brain imaging at single-cell resolution. She received her PhD from Tulane University in 2020. During her PhD study, her research focused on the development and optimization of light sheet fluorescence microscopy and its application in 3-D histopathology.
Guang Li is a PhD candidate in the Biomedical Engineering Department of Tulane University. He received his BS degree from Dalian University of Technology in China and his MS degree from Tulane. Working in the Translational Biphotonics Laboratory directed by Dr. J. Q. Brown, he is focusing on image processing of dual-view inverted selective plane illumination microscopy (diSPIM) and improvement of tissue clearing and fluorescent staining protocols. He also aims to apply light sheet microscopy to other practical areas.
Sharon E. Fox, MD, PhD, is the associate director of research and development in the Department of Pathology at LSU Health Sciences Center, New Orleans, and an anatomic pathologist at the Southeast Louisiana Veterans Healthcare System. She received her MD degree from Harvard Medical School in 2008 and her PhD from MIT in 2012. Her research interests include digital pathology, autopsy pathology, and the interaction between healthcare providers and image data.
Shadia I. Jalal, MD, is an associate professor of medicine in the Department of Medicine, Division of Hematology Oncology at Indiana University School of Medicine. She is a thoracic oncologist and phase I physician. She is a clinical researcher that designs, runs, and serves as the principal investigator of many clinical trials. Her research focus is in lung and esophageal cancers.
John Turek received his PhD from the University of Illinois at Urbana in veterinary medical science. He is a professor in the Department of Basic Medical Sciences at Purdue University College of Veterinary Medicine. For the past 16 years, his research has been in the area of biodynamic imaging with a focus on the use of 3-D cell cultures and BDI for derisking drug discovery and for cancer chemotherapy selection.
J. Quincy Brown is the Paul H. and Donna D. Flower Assistant Professor in the Department of Biomedical Engineering at Tulane University, and a program member in the Tulane Cancer Center. His research interests include point-of-care digital pathology for cancer diagnosis. He is a member of SPIE, OSA, IEEE-EMBS, and BMES.
David D. Nolte is the E.M. Purcell Distinguished Professor of Physics at Purdue University. He received his BS degree from Cornell University in 1981 and his PhD from the University of California at Berkeley in 1988, followed by a postdoctoral appointment at AT&T Bell Labs. He has been elected as a fellow of the Optical Society of America, a fellow of the American Physical Society, and a fellow of the AAAS. He is the technical founder of two biotech startup companies.