1.

Introduction

Tumor heterogeneity presents a challenge for the successful treatment of cancer using chemotherapeutics.¹ 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 microenvironment^5,6 and epigenetic variations^5,7–9 pose significant challenges to the selection of treatment based on genetic profiles. This has led, as an alternative, to phenotypic profiling^10–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 holography^15–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)²⁰ 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)²¹ 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.²² 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.²⁵

The methodology of TDS is usually applied to entire samples that can be as large as $1{\mathrm{mm}}^3$ 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.

Fig. 1

(a) An OCI image of a human esophageal biopsy and differential spectrograms (defined below) for the two circled areas. At time $t=0$ (black line), 100 μL of $25\mu \mathrm{M}$ cisplatin and $25\mu \mathrm{M}$ fluorouracil (5fu) combined solution is added to the sample. The two regions have significantly different responses: region 1 shows enhancement in low frequency and suppression in high frequency (redshift), while region 2 has suppression across all frequencies. (b) The terminal spectra of regions 1 and 2, respectively, compared to the average sample baseline spectrum.

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,²⁶ full-field OCT,²⁷ or conventional OCT^28–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.

2.

Materials and Methods

2.2.

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 $\sim 10\mu \mathrm{m}$. A QImaging EMC2 camera is used for image acquisition. Lenses L5 and L6 constitute a $4f$ 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 $256\times 256$ array.²¹ Because of the use of long focal lengths ($f_5=f_6=150\mathrm{mm}$), the speckle size on the camera plane is $\sim 60\mu \mathrm{m}$, and the reconstructed object-plane point-spread function is $\sim 15\mu \mathrm{m}$ spanned by three pixels.

Fig. 2

A schematic of the interferometer with a Mach–Zehnder off-axis digital holography configuration. The camera is on the Fourier plane. Optical sections are reconstructed using a 2-D spatial Fourier transform. ND, neutral density filter; ${\mathrm{L}}_1-{\mathrm{L}}_7$, lenses; ${\mathrm{BS}}_1-{\mathrm{BS}}_3$, beam splitters; ${\mathrm{M}}_1-{\mathrm{M}}_5$, mirrors; FP, Fourier plane; IP, image plane; ${\mathrm{f}}_5={\mathrm{f}}_6=150\mathrm{mm}$; ${\mathrm{f}}_7=50\mathrm{mm}$. The translation stage defines the coherence gate for time-domain ranging.

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.³⁴ 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 $(i,j)$ in the sample is calculated as the square of the FFT of the intensity time series:

Eq. (1)

$$S(i,j;f)={|{\mathcal{F}}_t[I(i,j;t)]|}^2.$$

The average spectrum of a region $\sigma$ is calculated as

Eq. (2)

$$S_{\sigma ,\text{raw}}(f)={⟨S(i,j;f)⟩}_{i,j}=\frac{\sum _{(i,j)\in \sigma }S(i,j;f)}{\sum _{(i,j)\in \sigma }1},$$

where $\sigma$ is a spatial mask segmenting the entire sample. The “raw” spectrum is normalized based on Parseval’s theorem:³⁴

Eq. (3)

$$S_{\sigma ,\text{norm}}(f)=\frac{S_{\sigma ,\text{raw}}(f)}{\sum _f{S}_{\sigma ,\text{raw}}(f)},$$

i.e., dividing the raw spectrum by the sum of all frequency components (including the DC component). The full-length spectrogram of segment $\sigma$ is a time-lapse sequence of differential spectra where each spectrum at time $\tau$ is calculated as a normalized logged spectrum subtracted by the baseline spectrum averaged over the full sample:

Eq. (4)

$$\mathrm{d}{S}_{\sigma }(f;\tau )=\log S_{\sigma ,\text{norm}}(f;\tau )-\log \frac{1}{N}\sum _{i=1}^N{S}_{\sigma ,\text{norm}}(f;\tau ={\tau }_i),$$

where $N$ is the number of baseline loops.

2.3.

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).³² 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,²² and the features of a spectrogram are represented by a vector of feature values.

Fig. 3

(a) A subset of the spectrogram masks used in the color merge. (b) Two maps of drug response of an esophageal sample “150903-14” (same as in Fig. 1) exposed to cisplatin and fluorouracil combination therapy under masks G0 and G1. (c) A “merged” bivariate color image with its 2-D color map.

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 as

Eq. (5)

$$\mathrm{d}{S}_{\sigma (\tilde{i},\tilde{j})}(f;\tau )=\log S_{\sigma (\tilde{i},\tilde{j}),\text{norm}}(f;\tau )-\log \frac{1}{N}\sum _{i=1}^N{S}_{\mathrm{\Sigma },\text{norm}}(f;\tau ={\tau }_i),$$

where the average baseline spectrum of the entire sample is used instead of the TDSI pixel’s “own” baseline. For a given mask, a feature value (discussed above) can be calculated for each pixel, and the feature values of an entire sample produce a 2-D image called a TDSI. Two TDS images under the “G0” mask and “G1” mask (blueshift associated with increased intracellular speeds) are shown in Fig. 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 cartography^36,37 and medical imaging,³⁸ 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 $8\times 8\text{pixel}$ averaging (instead of the standard $2\times 2\mathrm{px}$). The coarse eight-pixel averaging is used only to calculate the heterogeneity biomarker values of a sample at a scale of $40\mu \mathrm{m}$ and above. Finer scale is not necessarily meaningful for quantifying spatial heterogeneity across a millimeter sample. The feature values are bounded to a range $[-A_{\mathrm{th}},A_{\mathrm{th}}]$ and then assigned “scores” ranging from 0 to 1 for both heterogeneity evaluations, calculated using the following steps:

• 1. Select a set of $n$ masks

• 2. For each mask $u$, calculate $\mathrm{\Delta }{a}^u$ and $\mathrm{\Delta }[\mathrm{sgn}(a^u)]$

• 3. For each mask pair $u\text{-}v$, calculate $|\rho (a^u,a^v)|$ and $|\rho [\mathrm{sgn}(a^u),\mathrm{sgn}(a^v)]|$

• 4. The first biomarker denoting overall intramask heterogeneity is calculated as

  Eq. (6)

  $$h_1=\frac{1}{2n}[m_1\sum _{i=1}^n\mathrm{\Delta }{a}^u+m_2\sum _{i=1}^n\mathrm{\Delta }\mathrm{sgn}(a^u)].$$

• 5. And the second biomarker representing overall intermask heterogeneity is calculated as

  Eq. (7)

  $$h_2=1-\frac{1}{n(n-1)}\sum _{i-j\text{pairs}}\{|\rho (a^u,a^v)|+|\rho [\mathrm{sgn}(a^u),\mathrm{sgn}(a^v)]|\},$$
  where $a^u=\{a^u(\tilde{i},\tilde{j})\}$ are values of the TDS image of mask $u$, $\mathrm{sgn}(a^u)$ is a map of signs of $a^u$, $\mathrm{\Delta }{a}^u$ is the standard deviation of $a^u$, $\rho (a^u,a^v)$ is the correlation coefficient of $a^u$ and $a^v$, and $m_1$ and $m_2$ are normalization factors

  Eq. (8)

  $$\{\begin{array}{l}{m}_1=\frac{1}{A_{th}}\\ m_2=1\end{array},$$
  based on Popoviciu’s inequality on variances.⁴³ We use $A_{\mathrm{th}}=0.3$ and $n=3$, 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: $h_1=0$ when $n$ TDS images are completely uniform (totally homogeneous), and $h_1=1$ 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 $h_2=0$. In the opposite case when all values are completely correlated then $h_2=1$. Examples will be provided in the next section to illustrate these two heterogeneity benchmarks.

3.

TDSI Results

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.

Fig. 4

(a) Bivariate color representation of drug responses of two samples treated with different drugs, showing two univariate maps and a “merged” bivariate color map. The first sample was refreshed with DMSO, while the second sample was treated with $25\mu \mathrm{M}$ fluorouracil (5fu). (b) Global and regional spectrograms of sample “151208-6” from (a). The global spectrogram has a relatively weak response ($\max <10\%$), while the two circled areas have 30% to 60% enhancement or suppression. Drugs were added at $t=0$ (black line on spectrograms).

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.

Fig. 5

More examples of bivariate TDS images showing sample-to-sample variability in drug responses. Masks are designated in the lower left corners of images, while lower right corners designate drug treatments. Drug abbreviations: DMSO, 0.1% DMSO in growth medium (used as a negative control); cisp, $25\mu \mathrm{M}$ cisplatin; 5fu: $25\mu \mathrm{M}$ fluorouracil; tax, $5\mu \mathrm{M}$ paclitaxel; carbo, $25\mu \mathrm{M}$ carboplatin. “+” indicates a combination of two drugs. The color map and scale are the same as in Fig. 3.

There are roughly three types of heterogeneity, shown in Fig. 6 along with their $h_1$ and $h_2$ 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.

Fig. 6

Three types of samples that have different levels of same-mask heterogeneity and cross-mask heterogeneity, with scores on the right.

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.

Fig. 7

Time-lapse TDS image of samples responding to drugs. Sample “170317-9”: response of a DLD spheroid sample treated with $10\mu \mathrm{M}$ nocodazole [same as in Fig. 3(d)] with the G2 mask, showing a silent core shortly after drug was added (0.7 h), which was “invaded” by the drug and later achieved a spatially homogeneous response (2.7 to 8.7 h). Sample “170606-15”: response of an esophageal biopsy sample in the control medium, also under the G2 mask.

4.

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.

Fig. 8

TDSI map, iSPIM 3-D-reconstructed image and H&E histology image for an esophageal biopsy sample. All scale bars are $100\mu \mathrm{m}$. (a) Bivariate colormap representing G0 (broadband inhibition) and G1 (blue shift). (b) Bivariate TDSI map of human esophageal biopsy responding to 5-fluorouracil (5-fu) with G0 and G2 as the two masks. The color map and color scale are the same as in Fig. 3. (c) The iSPIM image with DRAQ 5 (blue) and Eosin (pink) for the same sample in (a), and (d) histology image. The orientations of the images are not registered.

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 ($\sim 15\mu \mathrm{m}\text{voxel}$ 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.

5.

Discussion

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 $h_1$ and $h_2$ 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.²⁶ 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.