Inverse-power-law behavior of cellular motility reveals stromal– epithelial cell interactions in 3D co-culture by OCT fluctuation spectroscopy

The progression of breast cancer is known to be affected by stromal cells within the local microenvironment. Here we study the effect of stromal fibroblasts on the in-place motions (motility) of mammary epithelial cells within organoids in 3D co-culture, inferred from the speckle fluctuation spectrum using optical coherence tomography (OCT). In contrast to Brownian motion, mammary cell motions exhibit an inverse power-law fluctuation spectrum. We introduce two complementary metrics for quantifying fluctuation spectra: the power-law exponent and a novel definition of the motility amplitude, both of which are signal- and position-independent. We find that the power-law exponent and motility amplitude are positively ( p <0.001) and negatively ( p <0.01) correlated with the density of stromal cells in 3D co-culture, respectively. We also show how the hyperspectral data can be visualized using these metrics to observe heterogeneity within organoids. This constitutes a simple and powerful tool for detecting and imaging cellular functional changes with OCT.


INTRODUCTION
A richness of functional information can be garnered from the time-dependent properties of tissues revealed by speckle fluctuations in coherence imaging techniques, including intracellular motility [1,2], hemodynamics [3], apoptosis [4], and ciliary activity [5]. While speckle fluctuation models have been well-developed for particle flow and diffusion in OCT [6,7], we are only beginning to understand the fluctuation spectra arising from intracellular motions in living tissues. Tissue fluctuation measurements are a growing area of research for assessing drug-response in cancer models in vitro [8,9], oocyte viability [10], and cortical cell recovery after stroke in vivo [3].
Optical coherence tomography (OCT) is particularly well-suited for imaging optically thick 3D tissue culture models, which provide a microenvironment that recapitulates many in vivo features [11]. We have previously shown how OCT of 3D mammary epithelial cell (MEC) organoids provides the ability to non-invasively and longitudinally monitor the evolution of organoids over weeks [12]. We have also shown how coupling temporal fluctuation analysis with polarization-sensitive OCT in these models allows for preferential image contrast of MEC organoids against a background of diffusive gold nanorod probes [13]. In that study, we reported on the observation of a motility signal from MEC organoids with a decay constant of several seconds. We note that, in contrast to cellular migration studies, this apparent "motility" represents the rapid, in-place motions of intracellular components measured over a relatively short time scale (<5 min). In this paper, we report on a comprehensive study of this motility signal from MEC organoids without the use of any contrast agents. As part of this study we developed methods to scale up to large numbers of organoids to achieve sufficient statistics, including methods for unbiased segmentation of MEC organoids, stable parameters for describing the fluctuation spectra, and an intuitive visualization method for observing cellular-level heterogeneity within the organoids.
Here we apply these techniques to a stromal-epithelial cell signaling model relevant to breast cancer. Increasing evidence suggests that fibroblasts in the mammary stroma influence the development and progression of cancer in MEC [11,14]. MEC with a basal-like subtype (as defined by gene expression), such as those in this study, may have unique susceptibility to fibroblast signals that are both soluble and mechanical in nature [15]. It is noted that the basal-like subtype is a type of breast cancer associated with a particularly poor prognosis [16]. We have recently shown that the morphology of basal-like MEC organoids changes with progression (normal versus pre-malignant) and depends upon the presence of fibroblasts in 3D co-culture [12]. Here we study the fluctuation spectra of MEC in the same co-culture model to understand how the apparent motility is affected by stromal cell signaling.
As will be shown, our fluctuation spectra suggest that intracellular motility in this model is well-described by an inverse power-law, and is distinct from the Lorentzian distributions expected for diffusion. Interestingly, the inverse power-law (also known as 1/f noise) is scale-invariant and arises in systems with long-range correlations or memory [17]. The power-law exponent is well-behaved in these 3D cultures and provides new insight into the effect of fibroblasts on MEC motility.

A. 3D co-culture preparation
We used a 3D cell co-culture method previously described by Johnson et al [18], with specific culture methods described in our previous publication [12]. Briefly, MEC form a multicellular, polarized structure (an "organoid") in 3D culture that exhibits many features similar to acini in the in vivo human mammary gland [19]. In this study we employed two MEC lines: MCF10A ("normal"), and MCF10DCIS.com ("pre-malignant"), the latter of which form organoids that are similar to ductal carcinoma in situ (DCIS). MEC were cocultured with stromal cells consisting of hTERT-immortalized fibroblasts from reduction mammoplasty (RMF). All cells were maintained in 2D or 3D culture at 37°C and 5% CO 2 with culture medium consisting of DMEM/F-12 and added ingredients listed in [20]. 2D cultures were used to grow a sufficient cell population to seed 3D cultures. For 3D cultures, a 1:1 mixture of collagen I (2mg/mL):Matrigel (BD Biosciences) extracellular scaffold was prepared. 10 mm diameter microwells were pre-coated with 85 μL of scaffold before the addition of 180 μL cell-scaffold mixtures on top. After gelation, 250 μL of cell medium was added and changed every 2 days.
A variety of 3D co-culture conditions were prepared. Cultures with either normal MEC (n c =10, n o =46) or premalignant MEC (n c =10, n o =39) were prepared, where n c represents the number of independent cultures, and n o represents the total number of organoids imaged with that condition. While the seed density of MEC were kept constant at 30,000/cm 3 , the seed density of RMF were prepared to be 0 (n c =6, n o =25), 30,000/cm 3 (n c =6, n o =22), 90,000/cm 3 (n c =4, n o =21), and 300,000/cm 3 (n c =4, n o =17). Culture times from 8-9 days (n c =8, n o =23) and 14-18 days (n c =12, n o =62) were explored. After 19 days, some cultures were formalin fixed (n c =4, n o =27).

B. OCT system and data acquisition
We employed a spectral domain, polarization-sensitive OCT system, described in detail previously [21]. Briefly, the OCT system light source was a broadband Ti:Sapphire laser centered at 800 nm (Griffin, KMLabs, Inc.) providing a coherence length of 2.6 μm in air. A 30 mm focal length lens in the sample arm provided a transverse resolution of ~10 μm in aqueous media. For all studies, 6.0-7.5 mW of linearly polarized light (H) was incident upon the sample, and the co-polarized backscattered light (H) was combined with the reference beam and directed to a custom spectrometer. A 4096 pixel line scan camera (Dalsa Piranha) was used to capture spectral interferograms at an A-line rate of 2 kHz into the first 2048 pixels. Sets of 300 successive B-mode images were collected over 1.5 × 1.5 mm into 1000 × 1024 pixels in x × z, respectively, at a frame rate of 0.876 ± 0.004 Hz. We expect this frame rate to be sufficient to capture fluctuation dynamics of MEC organoids, which we have previously found to exhibit an autocorrelation decay time of 4.8 ± 2 s using the same OCT system [13]. Reference subtraction and a digital dispersion compensation method described previously [22] were employed to compute OCT images, where image pixel intensities I(x, z) were computed from the absolute value of the complex analytic signal obtained from the spectral-domain OCT data. Figure 1 displays examples from an OCT image stack collected of a 3D co-culture.

C. Image segmentation
Both manual and semi-automated segmentation methods were used to determine regions of interest (ROIs) for OCT fluctuation spectroscopic analysis. All image stacks were first divided into 3 groups of N=100 images (1-100, 101-200, 201-300 within the original stack) for independent analysis, in order to perform cross-validation. Manual segmentation was then used for the studies reported in Section 3.A below to establish the spectral processing metrics. For all subsequent studies, a semi-automated method was employed to provide unbiased identification of MEC organoids within each OCT image stack.
An overview of the semi-automated segmentation method is displayed in Figure 2. First, the stack-averaged OCT image was computed (Fig. 2a), and a pre-corrected image was generated using predetermined, row-dependent pixel offset and pixel signal amplitude curves. This pre-corrected image was displayed via a custom MATLAB GUI enabling manual selection of a region of ECM that excluded MEC organoids (Fig. 2b, red outline). By assuming a priori that the ECM comprises a region of spatially invariant contrast, the OCT image was then further corrected in depth to auto-balance the image contrast (Fig. 2c). Finally, a pixel intensity threshold of 2× the ECM level and a series of Gaussian blur and erosion operations were performed to identify potential ROIs; only ROIs greater than 800 pixels (corresponding to ~30 cells in cross-section) were kept and identified as MEC organoids (Fig. 2d, yellow regions). Importantly, other than the user identification of the ECM ROI, the method is fully automated. Only one ROI out of 86 was manually rejected due to an image artifact; otherwise all auto-segmented ROIs were included in the subsequent analyses.
Control 3D cultures, seeded with RMF only, were also imaged. The RMF were not easily detectible in the OCT images. Therefore, we assume that segmentation of MEC:RMF cocultures predominantly selects MEC organoids.

D. Speckle fluctuation spectroscopic analysis
The speckle fluctuation power spectral density at each image pixel I(x, z) was computed from the discrete Fourier transform, F, over time, t, according to S(x, z, f) = |F(I(x, z, t))| 2 . This hyperspectral data was then spatially averaged over each ROI, and the f ≤ 0 terms were omitted before model fitting. Nonlinear least-squares fitting was employed with weights given by 1/S(f) (assuming that the variance is equal to the power spectral density). Two models were fitted: a Lorentzian model, which is known to represent Brownian (diffusive) motions in OCT data [23], and an inverse power-law model, which has been used to describe self-organized critical systems [24]. The Lorentzian can be written as: (1) where the free-fit parameters were taken to be the damping term, γ, a constant of proportionality, c 0 , and additive white noise, n, all of which were constrained to be >0 during the fitting procedure. The inverse power-law model can be written as: (2) where the free-fit parameters were taken to be the exponent, α, a constant of proportionality, c 0 , and additive white noise, n, all of which were similarly constrained to be >0 during the fitting procedure.

E. Speckle fluctuation amplitude ("motility amplitude") analysis
In addition to the quantification of the spectral shapes via γ or α in the models above, a measure of the overall fluctuation amplitude provides a complementary metric for cellular function. While the proportionality constant c 0 in Eq. (1) and (2) might appear useful in this regard, it was found to vary strongly with the overall pixel intensity and thus is not representative of cellular function. In our previous work, we introduced the motility amplitude, M, as a normalized standard deviation as follows [13]: (3) where N is the number of images, t i are the discrete times at which each image in the stack is collected, and the overbar indicates a time average. Unfortunately, as illustrated in Fig. 3b, this metric does not provide a consistent value for organoids within the same culture condition. Despite normalization for the presumed effect of shot noise in the denominator of Eq. (3), M exhibits a strong decay with the depth of the organoid within the OCT image.
To address this shortcoming, here we propose a new method for estimating the motility amplitude based upon a modified standard deviation as follows: (4) where Γ I is the temporal autocorrelation of I(x, z, t), and Δt is the image sampling time.
Thus, we write . The autocorrelation at Δt quantifies the OCT signal that persists for one sampling time, which has the important effect of omitting random noise (such as shot noise) that decorrelates instantaneously. Thus, M represents the fluctuation amplitude of signals that decorrelate at times >Δt but shorter than the total measurement time. Because cell-driven fluctuations have been shown on average to decorrelate with a decay time of ~4.8 s [13] equivalent to ~6Δt, M effectively captures the cellular fluctuation signal.
For M to be unitless and independent of pixel intensity, normalization by the average pixel intensity is used in Eq. (4). Thus, M represents a fractional fluctuation amplitude. As shown in Fig. 3c, this new metric is extremely consistent across a wide dynamic range of organoid average pixel intensities within the same culture, and thus organoids physically lower in the image (which exhibit SNR roll off) provide an M consistent with those nearer the top of the OCT image. Quantitatively, this can be stated by noting that M was independent of both average pixel intensity and depth in the image (R 2 < 0.03 for both for the ten ROIs of Fig.  3), while intensity and depth were yet strongly inversely correlated (R 2 > 0.9).
For completeness, we also considered renormalizing Eq. (3), the standard deviation, similarly to Eq. (4) so that it is a unitless number. Specifically, the numerator of Eq. (3) was divided by the denominator of Eq. (4). The result of this test is displayed in Fig. 3d, which clearly demonstrates that it is inadequate to provide the depth-independent results needed to represent organoid motility (and, in fact, overcompensates the data, making the lower organoids now display a higher motility value). Thus, we conclude that the modified standard deviation proposed in Eq. (4) and shown in Fig. 3c represents a significant new metric for quantifying the amplitude of cell-driven speckle fluctuations. Importantly, regions of high motility amplitude are well correlated with the ROIs determined semi-automated segmentation, as displayed in Fig. S1. This definition of M is subsequently used in all of the studies described below, where M(x, z) is spatially averaged over each ROI.

RESULTS AND DISCUSSION
A. Evidence for the appropriateness of the power-law exponent (α) and the motility amplitude (M) as metrics for quantifying cellular function 1. α uniquely describes live cells-The existence of a reduced set of parameters that accurately represent OCT fluctuation spectra would constitute a powerful tool for 3D cell culture analysis. We first sought to determine whether the power spectra from MEC organoids were consistent with either the Lorentzian power spectrum (Eq. (1)) or the inverse power-law model (Eq. (2)). To compare the power spectra of living and non-living samples, we evaluated ROIs to include either ECM or MEC organoids; 3 of each were manually segmented from each image stack. ROIs were of approximately equal area, but were chosen to be of varying depth within the image to vary the OCT signal intensity. All cell culture conditions described in section 2.A. above were tabulated. The results of this analysis grouped by ROI type and fixation state are summarized in Fig. 4. Figure 4 shows the coefficient of determination (R 2 ) values obtained from the model fittings of fluctuation spectra within ROIs comprised of ECM or MEC organoids. Regions of ECM exhibited lower R 2 overall, and thus were less consistent with both models in comparison to organoids. This may be due to the fact that ECM has lower OCT signal intensity, causing its diffusive motions to lie below the noise level. Within organoids, both models exhibited R 2 distributions entirely above 0.90. Importantly, only the fluctuation spectra from live (unfixed) organoids exhibited a preference for the inverse power law model, with a median value of R 2 =0.99. The fact that formalin-fixed cells did not exhibit this preference suggests that the inverse power-law model only describes live cells. This is consistent with the picture that MEC exhibit self-organized criticality, i.e., tend toward a minimally stable state where small perturbations cause a scale-invariant response [24]. This scale-invariance, or fractal nature, can be seen by the fact that α in Eq. (2) is a unitless quantity. It is important to note that α is measured over a limited frequency band from 9-440 mHz in this study; in the absence of further data the inverse power-law description should not be extrapolated to all scales (a requirement for a true fractal object [25]). That said, in other areas of research, such as tissue mechanics, a fractal model is useful for describing viscoelasticity [26], which is typically employed in practice over a limited, experimentally-available frequency range [27].
2. Both α and M are independent of OCT signal and depth-One of the most challenging aspects of quantifying OCT images is the depth roll-off in SNR which arises from a combination of limited depth-of-focus, light attenuation, multiple light scattering [28], and in spectral-domain OCT, finite spectrometer resolution [29]. In addition, the OCT system signal strength from a given object can change from day-to-day. The robustness of the power-law exponent, α, and motility amplitude, M, were tested against these potential sources of systematic error as follows. To simulate the effect of the depth roll-off, the same MEC organoids were imaged as they were translated to five different depths in random order. To simulate the day-to-day variations in OCT signal, the reference arm power was attenuated to three different settings (nominally 100%, 75%, and 50% of maximum power) while imaging the same MEC organoids. These results were then compared against the endogenous variation in the organoids, i.e., the cross-validation results obtained when dividing image stacks into 3 groups of 100 images. For this study, two additional 3D cultures each were prepared with normal MEC or premalignant MEC, both without fibroblasts, and imaged at 1 week of culture time. The depth roll-off n o =6, and the OCT signal n o =7, where n o is the number of MEC organoids imaged.
The results of this study are summarized in Table 1 in terms of the standard deviation (as a percentage of the average value) of α and M each. Increased standard deviations due to the added variation (either by modifying the signal or the depth) are displayed in the bottom row. Initially, we found that the endogenous variation in α and M are ~3-4% and ~6-7%, respectively. When tracking the same MEC organoids at varying signal levels, the standard deviation did not significantly increase. However, varying the image depth did cause significant increases, particularly in the standard deviation of M to a value of 12.8%. Despite this effective doubling in the error of M from that endogenous to the organoids, in the context of the variability inherent to 3D tissue cultures, the metric M is yet sufficiently quantitative to evaluate cellular functional changes, as will be shown below.
3. Both α and M vary with fixation-To further validate that the proposed fluctuation spectra metrics represent cellular function, we measured α and M of MEC organoids before ("live") and after ("fixed") exposure to formalin fixative. The same 3D cultures as described in Section 3.A.2 were used, with n o =37 for both live and fixed.
The results of this study are summarized in Figure 5. Interestingly, α increases after cell fixation, which represents that the power spectral density is relatively lower at higher frequencies. It is important to remember that the inverse power law model fitting (R 2 values in Fig. 4) gets worse as cells are fixed; thus, the increase in α from ~1.6 to ~2.2 may represent more Lorentzian-like (Brownian) behavior, since a Lorentzian scales approximately as f −2 for small γ.
In contrast, the value for motility amplitude decreases after cell fixation, as expected. However, M does not exactly reach zero for fixed cells, where a non-zero offset of M≈0.11 represents the noise floor remaining after rejection of the instantaneously decorrelating noise via Eq. (4).
It is important to note that the fixation process stiffens collagenous tissue via cross-linking. This could cause some of the effects observed in Fig. 5 suppressing Brownian motions in the sample. Thus, the results of this experiment are necessary but not sufficient to show that α and M represent cellular function and not just Brownian motion. Taken together with the other evidence above (Fig. 4 showing the preference for the inverse power law only for live cells, and Fig. 3c showing excellent correspondence between regions with high motility and the positions of the organoids), it seems exceedingly likely that the signals represented by α and M represent primarily ATP-driven, and not Brownian, motions Furthermore, the use of speckle fluctuations for elucidating cellular responses to a variety of conditions has been well-established by other research groups [8,30,31]. As will be demonstrated below, these metrics are useful for indicating stromal-epithelial cell interactions in 3D co-cultures.

B. The power-law exponent reveals stromal-epithelial signaling
Interactions between stromal fibroblasts and mammary epithelial cells are implicated in breast cancer tumorigenesis, and 3D tissue co-cultures have been employed to model these interactions in controlled settings [32]. Here we hypothesized that α and M of MEC organoids would exhibit fibroblast-dependent changes in co-culture, similar to morphological changes that we previously observed in the same model system [12]. Using the large data set described in the methods (n o =85, 3 measurements each), we sought to determine the dependence of α and M on each of the variations explored, including the MEC cell line (normal vs. pre-malignant), the fibroblast seed density, and the culture time. We performed linear regression using R (v. 3.1.1) according to the following model: (5) where x RMF is the seed density of fibroblasts in 10 10 /m 3 , n MEC is a dummy variable defined as 0 for normal MEC and 1 for pre-malignant MEC, t is the culture time in days, m g is an indicator variable for group number within the stack (as an internal control), and β i are the free-fit parameters for which estimates are summarized in Table 2. Note that each measurement, including the 3 replicate measurements of each organoid, were treated as independent in this analysis.
The results of the linear regression underscore the importance of stromal fibroblasts in MEC motility. While the order of the repeat measurements (represented by the group number, β 4 ) did not have any significant impact on α or M, fibroblast seed density strongly predicted α and M (p=10 −14 for α, and p=10 −3 for M). Interestingly, increasing fibroblast density increases α (selectively suppressing higher frequency fluctuations) while decreasing M (suppressing the overall motility amplitude). Figure 6 illustrates this dependence.
The observed stromal-dependence of α and M may be an effect of increasing stiffness in the co-cultures from the fibroblasts, which was qualitatively palpable when working with the co-cultures. Importantly, we note that the co-culture gels were not allowed to detach from the container walls, and therefore shrinking of the gel volume cannot explain the increase in ECM stiffness. Rather, the fibroblasts increase tension and/or remodel the ECM thereby modulating stiffness. Mechanical changes in the ECM have previously been implicated in modifying MEC organoid behavior [33], and both soluble biochemical signals and mechanical signals from fibroblasts are known to independently affect MEC activity [15].
Interestingly, α was independent of the culture time and MEC line, while M exhibited a slightly significant (p<0.05) dependence on them (Table 2). Interestingly, pre-malignant cells exhibited higher motility amplitude than normal cells, and longer culture times tended to suppress M. The fact that the cell line is not strongly indicated is somewhat unexpected, given the strong cell line dependence of organoid morphology we previously observed in this co-culture system [12]. However, we note that in the previous study, we found that the morphological changes were larger only when both fibroblast seed density was larger and cell type was pre-malignant, suggesting an interaction between these variables. Also, our other previous work has suggested that the pre-malignant MEC line (MCF10DCIS.com) may be more susceptible to signals from fibroblasts in co-culture [32]. Therefore, in a second statistical analysis, we modeled the motility amplitude according to the following, which includes an interaction term for MECs and RMFs: (6) where the new fitting parameter, β 5 , represents the strength of the interaction between fibroblast density and MEC line. We omitted β 4 (group number) since it was previously not significant and to limit the number of parameters in the linear regression. The results are summarized in Table 3. Using a linear model with an interaction term, Eq. (6), we find that indeed, the interaction between fibroblast seed density and MEC line is significant (p=10 −3 ), with a higher density of fibroblasts and the pre-malignant cell line being correlated with a higher motility amplitude, M. This is interesting considering that the linear term for fibroblast seed density, β 1 , is negative which indicates an overall suppression of M with increasing fibroblast concentration. Thus, while the overall activity of fibroblasts is to suppress motility, premalignant MEC organoids are differentially modulated by their presence, exhibiting larger M in their presence than do normal MEC organoids. This finding is consistent with the previous morphological study [12], and thus we propose that the model of Eq. (6) is more appropriate than that of Eq. (5) to describe M in this co-culture system.

C. Visualization of parameterized hyperspectral images
One of the advantages of OCT is the ability to visualize the spatial heterogeneity of cellular function in 3D cultures. MEC organoids are known to develop into polarized, acinar structures over time which would be expected to cause heterogeneity in the observed power spectra. Importantly, the effect of fibroblasts from discrete locations around organoids may cause heterogeneous changes in ECM structure or tension to modify the behavior of individual cells within the organoids. Current genomic profiling methods for elucidating stromal-epithelial interactions cannot readily separate genes from epithelial and stromal cells in co-cultures (or in real tissues), and cannot reveal spatial heterogeneities within organoids. Thus, visualization of OCT-based motility metrics might provide insight that is complementary to gene expression profiling.
In traditional hyperspectral imaging [34] or in spectroscopic OCT [35], individual parameters of interest (such as wavelength, or tissue oxygenation) are often displayed on the original 2D space via a straightforward color map. Here we have reduced the temporal speckle fluctuation spectra to two parameters, α and M, to describe the spectral shape and amplitude, respectively. Because the ability to describe the spectral shape by α depends on the goodness-of-fit to the inverse power-law model (Eq. (2)), the R 2 of the fitting is a third parameter that should also be monitored.
When there is a need to display multiple scalar fields (such as α, R 2 , and M) on the same 2D surface, sparse-glyph techniques are often preferred [36]. Here we employed the opensource platform ParaView to render these 3 motility metrics onto the original B-mode OCT image space. Because spatial averaging is needed to suppress the effects of noise on α, power spectra were averaged over a pillbox region of 9 pixels in diameter centered at each pixel within the ROIs (where pixels outside the ROI were omitted from the pillbox). This pillbox diameter is equivalent to 13.5 μm, or a cross-sectional area equivalent to that of ~2.5 MEC [12]; α and the associated R 2 were then computed for each spatially-averaged spectrum, while M was computed for all pixels without spatial averaging. M was then displayed as a grayscale image, with contours from the ROIs drawn in yellow. Similar to a technique presented in [37], α was presented as a glyph of varying color. The color was linearly mapped to a double-ended blue-red color map; blue being chosen to represent small values of α that have relatively higher (blue-shifted) frequencies, and red being chosen to represent larger α that have relatively lower (red-shifted) frequencies. Importantly, the sizes of the glyphs were scaled by the R 2 value to indicate greater confidence as larger glyphs. The scaling was computed according to the following nonlinear expression: (7) where d is the fractional diameter of the glyph, and R 0 2 is the target value of R 2 at a fractional diameter of 50%. Here we set R 0 2 = 0.8 to emphasize that R 2 <0.8 indicates a poor model fit.
The results of this visualization method are illustrated in Figure 7, where we observe that there are wide variations in motility over the MEC organoids that are not captured by the ROI-averaged analyses above. First, we note that the ROIs themselves (yellow contours) do not always represent regions of high M (white background). Although the ROIs represent regions of high light scattering from the time-averaged OCT image (Fig. 2d), clearly there are regions with high OCT signal that exhibit low motility amplitude, which could be inactive cells within the organoid. Conversely, there are also regions with low OCT signal that exhibit high M, which could represent regions that the MEC are actively extending into. We also find that there is an apparent correlation between organoid regions of low M and red-shifted power spectra (high α), an example of which is apparent in Fig. 7c, suggesting that the possible inactive cells exhibit both low M and red-shifted power spectra. Interestingly, regions of high M within organoids exhibit values of α ranging from red-to blue-shifted with varying levels of confidence (R 2 ). Because the α values represent the response of <3 cells each, it is possible that isolated red-shifted glyphs with high M represent a small cluster of cells with a different functional state.

CONCLUSIONS
Here we present biophysical evidence of stromal-epithelial interactions in 3D breast tissue co-cultures via monitoring speckle fluctuations within the 9-440 mHz band. Significant changes in the α and M of MEC organoids are observed as a function of fibroblast seed density. The signal-and position-invariance of these metrics provide a powerful tool for measuring and visualizing cell function in 3D cultures; when using OCT systems with the same wavelength and axial resolution, we expect that these metrics will be consistent across laboratories. Coupling these techniques with OCT-based morphological analysis of the MEC organoids [12] provides a comprehensive picture of stromal-epithelial signaling. In particular, we found a similar response of the motility amplitude, M, to interactions between fibroblasts and the MEC line (normal vs. pre-malignant) as we previously did with morphological metrics including organoid size and asphericity. The ability to visualize the spatial heterogeneity of the motility of small numbers of MEC (in this study, <3 MEC per averaging window) may lead to new insight such as locating cancer stem cells [38] that cannot be resolved from gene analysis. The intuitive visualization method we propose here may allow one to utilize the already effective human visual system to identify spatial trends in organoid behavior. Importantly, all of these studies can be accomplished non-invasively, without exogenous labels, and longitudinally, on live 3D cultures.
An open question remains as to the underlying mechanism of the fluctuation signals, which fundamentally arise from non-Brownian motions of light scattering components of the cells. The relevance of the inverse power-law in this study implicates biophysical processes with scale invariance (i.e., self-similarity as described by fractals) and long-range memory, although the extension of these measurements to at least 3 decades in frequency would be needed to unambiguously prove this hypothesis [25]. The particular cellular components contributing to the light scattering signal can vary widely depending upon the cell type, where organelles such as the nucleus and mitochondria are often implicated [39]. We also expect that possible changes in cellular morphology may underlie some of the observed changes in the dynamic light scattering signals. Bulk cellular velocity may also be expected to play a role, as linear velocities of MEC in 3D culture on the order of 0.2μm/min [40] would be expected to provide speckle fluctuations at the low end (~20 mHz) of our measurement band. There is growing interest in characterizing the non-equilibrium and molecular motor-driven motions of subcellular components [41][42][43] which, in future work, could be used to model and better understand the observed OCT fluctuation spectra. Despite these unanswered questions, efforts to profile speckle fluctuations phenotypically can provide insight and predictive value for applications such as drug screening [8], and our efforts here demonstrate its utility for stromal-epithelial systems.
Importantly, stromal-epithelial signaling affects the progression of basal-like breast cancers, a subtype of breast cancers associated with a poor prognosis. The methods developed here represent a comprehensive tool for understanding the interactions of these cells in the challenging imaging environment of 3D co-culture models that recapitulate the in vivo microenvironment.

Supplementary Material
Refer to Web version on PubMed Central for supplementary material.  Example of the semi-automated segmentation method (see text for description) applied to a stack-averaged OCT image of MEC organoids (pre-malignant) in 3D co-culture with RMF (seed density 90,000/cm 3 ) at 14 days culture time. Notably, these MEC organoids are significantly less spherical than those shown in Fig. 1 due to the different cell line (premalignant vs. normal), increased RMF seed density, and increased culture time, consistent with our previous findings on MEC organoid morphology [12]. The challenging MEC organoid shapes are reliably segmented into 10 ROIs in this example using the semiautomated method.  Box-and-whisker plots of the distribution of R 2 values obtained from least-squares fitting to the Lorentzian and inverse power-law models for all cell culture conditions. Only live cells exhibited a preference for the inverse power law model (red arrow). (ECM unfixed n R =180, ECM fixed n R =36, MEC unfixed n R =153, MEC fixed n R =33, where n R is the number of manually-selected ROIs of each condition).  Box-and-whisker plots of the distribution of α and M at each fibroblast seed density (for both cell lines and all culture times). As seed density increases, α increases strongly, while M decreases slightly. The bottom panel illustrates inverse power-law model fittings of the power spectra for representative ROIs seeded with no fibroblasts (blue) compared to 300 μL −1 fibroblasts (purple); the noise constant (n in Eq. (2)) was subtracted and spectra normalized before plotting to emphasize the effects of α. As shown, a higher value of α indicates the suppression of higher frequencies relative to lower frequencies. Visualization of motility metrics on the same data as in Figs. 2 and 3, where panels (b) and (c) are close-up views of organoids (yellow arrows) in panel (a). The background grayscale image indicates M, and α is overlaid as a spherical glyph colored from blue to red to represent fluctuation spectra with more high frequency or low frequency components, respectively. Blue or red hollow glyphs indicate α values beyond the scale (<1 and >2, respectively). The size of the glyph represents R 2 from the power-law fitting, with larger glyphs indicating α values with greater confidence. Yellow contour lines indicate the ROIs determined from segmentation of the time-averaged OCT image.  Table 3 Parameter Estimates for Motility Amplitude According to Eq. (6)  Estimate (standard error) obtained using linear regression according to Eq. (6).