## Abstract

Altered mechanical properties of the tumor matrix have emerged as both the cause and consequence of breast carcinogenesis. Increased tumor stiffness has traditionally provided a viable metric to screen for malignancies via palpation or imaging. Previous studies have demonstrated that the microscale mechanical properties of the cell substrate influence tumor proliferation and invasive migration *in vitro*. Nevertheless, the association of the mechanical microenvironment with clinical hallmarks of aggressiveness in human breast tumors, including histopathological subtype, grade, receptor expression status, and lymph node involvement is poorly understood. This is largely due to the lack of tools for mapping tumor viscoelastic properties in clinical specimens with high spatial resolution over a large field of view (FoV). Here we introduce laser Speckle rHEologicAl micRoscopy (SHEAR) that for the first time enables mapping the magnitude viscoelastic or shear modulus, |G*(x,y,ω)|, over a range of frequencies (ω = 1–250 rad/second) in excised tumors within minutes with a spatial resolution of approximately 50 μm, over multiple cm^{2} FoV. Application of SHEAR in a cohort of 251 breast cancer specimens from 148 patients demonstrated that |G*(x,y,ω)| (ω = 2π rad/second) closely corresponds with histological features of the tumor, and that the spatial gradient of the shear modulus, |∇|G*(x,y,ω)||, is elevated at the tumor invasive front. Multivariate analyses established that the metrics, (|G* |) and (|∇|G* ||), measured by SHEAR are associated with prognosis. These findings implicate the viscoelastic properties of the tumor microenvironment in breast cancer prognosis and likely pave the path for identifying new modifiable targets for treatment.

Laser speckle rheological microscopy establishes the links between microscale heterogeneities of viscoelasticity and histopathological subtype, tumor grade, receptor expression, as well as lymph node status in breast carcinoma.

## Introduction

It is increasingly recognized that the pathogenesis and progression of breast carcinoma are guided by the mechanical transformation of the tumor microenvironment (1–3). *In vitro* studies demonstrate that altered mechanical properties of the extracellular matrix (ECM) promote neoplastic transformation, differentiation, growth, and migration of mammary cells (2–4). Although these insights are well established through *in vitro* mechanobiology studies, their implications for the clinical course of disease remain unclear (5–8). Breast carcinomas are miscellaneous, categorized by their histological and molecular subtypes, grade, and lymph node status (9, 10). These prognostic criteria determine aggressiveness, predict response to therapy, and suggest modifiable targets for treatment (10). The relationships between mechanical remodeling of tumors and these prognostic criteria are poorly understood. If the link between micromechanical features and clinical phenotypes is identified, prognostic classification of breast malignancies may be improved. Moreover, novel therapies could be developed on the basis of modifying the micromechanical drivers of carcinogenesis.

Increased stiffness or elasticity, a consequence of desmoplastic reaction, has been linked to metastasis and drug resistance in solid tumors (2, 11–15). However, breast tissue is not merely elastic, but rather exhibits both viscous and elastic behaviors under different loading conditions, best characterized by the frequency-dependent complex viscoelastic modulus also termed the shear modulus, G*(ω) (16). Emerging evidence suggests that tumor cells sense both the microscale elastic and viscous traits of their substrate and respond to viscous dissipation by either activating or suppressing the oncogenic reprogramming, in ways not explained by changes in elasticity alone (16, 17). Moreover, spatial heterogeneities of |G*(ω)| within the tumor may create gradients, that is, **|**∇|G*(ω)|**|**, that foster directional migration and invasion (7).

Histopathological analyses reveal that breast cancer subtypes exhibit distinct tumor and stromal composition and architecture (9, 10). Deposition and rearrangement of the ECM components by neoplastic and host cells in turn cause mechanical heterogeneities, lending credence to the premise that the tumor mechanical microenvironment is associated with the histological subtypes (1). Mechanical properties of the ECM may regulate mitotic activity and subsequently govern histological grade (11, 14, 18). Expression of estrogen, progesterone, and HER2 (ER/PR/HER2) receptors in tumor cells are associated with fibroblast and immune cell infiltration. The subsequent fibrotic and inflammatory responses in turn alter the mechanical microenvironment of the tumor and may link viscoelastic properties to molecular expressions (1, 19–21). Given that preferential invasion of tumor cells is orchestrated by perpendicularly aligned ECM fibers at the tumor–ECM interface, also termed the invasive front (IF), mechanical heterogeneities are expected to correlate with lymph node involvement (7, 22, 23). Despite growing evidence, the links between tumor viscoelasticity and these prognostic criteria are missing and the potential implications for diagnosis and treatment remain distant, largely due to absence of tools for mapping the |G*(ω)| with high resolution.

Here, we introduce laser Speckle rHEologicAl micRoscopy (SHEAR) that estimates and maps the shear modulus, G*(x,y,ω), of tumors with high spatial resolution over large fields of view (FoV). SHEAR is developed by modifying an inexpensive upright microscope to illuminate the specimen and capture time series of speckle patterns, formed by multiple scattering and interference of coherent light from tissue particles (Fig. 1A and B). Thermal motion of intrinsic light-scattering particles within tissue creates dynamic phase shifts and speckle intensity modulations at a rate commensurate with the extent of particle displacements and in turn the viscoelastic properties of tissue microenvironment (Fig. 1C–E). Consequently, in breast, softer adipose and glandular tissues exhibit rapid speckle intensity fluctuations compared with the fibrous ECM. Therefore, by measuring the mean square displacement, MSD, of scattering particles from temporal speckle intensity fluctuations, and estimating the radius of scattering particles, the frequency-dependent complex shear modulus, G*(ω), can be assessed over a range of angular frequencies, ω. SHEAR is similar to dynamic light scattering (DLS) in that it measures the G*(ω) from the thermal fluctuations of scattering particles (24). However, unlike DLS, SHEAR is not limited to single scattering events in transparent samples and instead exploits multiple scattering from an ensemble of intrinsic light scatterers to measure heterogenous, stiffer turbid samples including tissue.

We demonstrate that SHEAR can rapidly map the magnitude of the complex shear modulus, |G*(x,y,ω)| with 50 μm spatial resolution over scalable FoVs of multiple cm^{2}, covering tumor epithelium, IF, and peritumoral ECM within minutes. The |G*(x,y,ω)| maps plotted at a frequency of ω = 2π rad/second closely corresponded with histomorphological features and elevated gradient of shear modulus, **|**∇|G*(x,y,ω)|**|**, observed at the IF. Our study shows for the first time that shear modulus and the indices of mechanical heterogeneity estimated by SHEAR are associated with prognostic indicators of histological subtype, receptors' expression status, grade, and lymph node involvement. These results establish the utility of SHEAR technology for studying the viscoelastic behavior of breast cancer and likely open impactful avenues for integrating these new micromechanical metrics in the diagnosis, treatment planning, and development of novel therapeutics.

## Materials and Methods

### SHEAR setup

A schematic of the SHEAR optical instrument is displayed in Fig. 1A. A laser beam (633 nm, Helium-Neon, 45 mW, JDSU) was focused at the back focal plane of the objective (10× Olympus, NA = 0.25), creating a 1.6 mm diameter beam (5 mW) at the sample surface. Backscattered time-varying speckle patterns (Fig. 1B) were imaged via the objective and tube lens (*f* = 175 cm) by a CMOS camera (Basler, acA 2000-340 km) operated at 250 frames per seconds (fps). The sample was placed on a motorized stage interfaced with a step-motor controller (ESP 301-3G, Newport Corporation). Speckle frame series were acquired over a region of interest (RoI) of 500 × 500 μm and the sample was translated in 450 μm steps in a serpentine pattern, maintaining a 10% of overlap between acquired RoIs to scan the entire sample. For each RoI, speckle frames were acquired for 1 second at 250 fps. Brightfield images were acquired simultaneously to facilitate coregistration with histopathology. A Custom-built C++ software was used for instrument automation to synchronize the motion controller and camera trigger.

### SHEAR image reconstruction

Methods detailed in our prior studies were extended to estimate the frequency-dependent shear modulus, G*(x,y,ω), at each pixel location, (x,y), from time-varying speckle frames (25–30). The speckle intensity autocorrelation, g_{2}(x,y,t), was measured from the speckle time series and the MSD of local scattering particles was subsequently calculated (Supplementary Data S1; Supplementary Fig. S1A and S1B). The G*(x,y,ω) was determined by substituting the MSD in the generalized Stokes–Einstein relation (GSER) via |${{\rm{G}}^{\rm{*}}}( {{\rm{x}},{\rm{y}},{\rm{\omega }}} ) = \frac{{{{\rm{K}}_{\rm{B}}}{\rm{T}}}}{{{\rm{\pi }}a{\rm{i\omega }}{\cal F}\{ {\langle {\Delta {{\rm{r}}^2}( {{\rm{x}},{\rm{y}},{\rm{t}}} )} \rangle } \}}}$|, where |$\ {\cal F}$| is Fourier transform and *a* is the scattering particle radius (Supplementary Data S1; Supplementary Fig. S1C; refs. 25–29, 31, 32). The particle radius, *a*, was estimated by calculating the ratio of the decay rates of g_{2}(t) curves acquired at perpendicular and parallel polarizations (Supplementary Data S2; Supplementary Fig. S2A–S2D). The magnitude of the shear modulus, |G*(x,y,ω)|, was evaluated over a frequency range, 1<ω<250 rad/second, determined by the acquisition time (1 second) and the camera frame rate (250 fps). Although |G*(x,y,ω)| was calculated over the entire range, 1<ω<250 rad/second, for ease of comparison with histopathology and correspondence with rheometry, in this article we display |G*(x,y)| maps at a single frequency ω = 2π rad/second (*f* = 1 Hz). From the |G*(x,y)| maps, the following metrics were calculated: the spatially averaged |$\overline {| {{G^*}} |}$| over all pixels (*x,y*); the shear modulus gradient, **|**∇|G*(x,y)|**|**; the spatially averaged gradient, |$\overline {| \nabla |{G^*}||}$| ; heterogeneities in |$\overline {| {{G^*}} |}$| in the tumor epithelium versus stroma; and the information entropy, H, of the |G*(x,y)| maps (Supplementary Data S4; Supplementary Fig. S3). Rewriting the shear modulus as G*(ω) = G'(ω)+iG"(ω), further permitted extracting the elastic, G'(ω), and viscous, G"(ω), moduli as the real and imaginary components (Supplementary Data S1). Figure 1D displays representative g_{2}(t) and MSD (inset) for a normal fibroadipose breast tissue. The frequency-dependent |$\overline {| {{G^*}( {\rm{\omega }} )} |} ,\ \overline {G'( {\rm{\omega }} )}$|, and |$\overline {G''( {\rm{\omega }} )}$| obtained by SHEAR are shown in Fig. 1E: All three curves closely agree with the results of conventional rheometry at frequencies below ω = 20π (62.8) rad/second, (i.e., 10 Hz), above which, the accuracy of rheometer measurements is influenced by the inertial effects as discussed below.

### Study design

The study was approved by the Institutional Review Board of Massachusetts General Hospital (MGH IRB#2011P000301; Boston, MA). Fresh, excess, and deidentified human breast specimens (*N* = 251) were collected from 126 patients undergoing breast cancer surgery (mastectomy or lumpectomy, 1 male, 125 females) and 22 patients undergoing breast reduction surgery, between March 2015 and October 2019. The clinical diagnosis for each specimen, the tumor size, histopathological grade, hormonal receptor status, lymph node status, and history of neoadjuvant chemotherapy or radiotherapy was obtained from pathology records. Table 1 summarizes the tumor diagnoses. The samples were refrigerated at 4°C and imaged within 2 hours of being grossed. Prior to imaging, the sample was marked with ink, warmed to 37°C in a water bath, and loaded onto the motorized linear stage within the SHEAR device (Fig. 1A). After imaging, each sample was fixed in 10% formalin, sectioned at 100 μm increments in depth (7-μm–thick sections), and stained with hematoxylin and eosin (H&E) and picrosirius red (PSR). The H&E slides were digitized using Nano-zoomer whole slide scanner (2.0 HT). PSR slides were digitized using a circularly polarized light microscope (BX43, Nikon, 4×). Fiducial ink marks in the sample photographs and brightfield images were used to coregister the SHEAR data with histology. Pathologist (E.F. Brachtel) blinded to the SHEAR data provided histopathological analysis from H&E slides.

Diagnosis/type . | Number . | |${\overline {| {\rm {G^*}} |}}$| (kPa) . |
---|---|---|

Histological type | ||

Invasive ductal carcinoma (IDC) | 75 | 8.8 ± 8.5 |

Invasive lobular carcinoma (ILC) | 15 | 10.5 ± 5.3 |

Invasive ductal lobular carcinoma (IDLC) | 7 | 5.9 ± 3.7 |

Invasive mucinous carcinoma (IMC) | 3 | 4.4 ± 1.4 |

Invasive papillary carcinoma (IPC) | 8 | 2.4 ± 1.3 |

Invasive micropapillary carcinoma (IMPC) | 3 | 19.7 ± 12.6 |

Histological grade | ||

Grade 1 | 10 | 13.9 ± 12.2 |

Grade 2 | 46 | 10.2 ± 8.7 |

Grade 3 | 55 | 6.3 ± 5.3 |

Estrogen receptor (ER) status | ||

ER^{+} | 101 | 9.2 ± 8.1 |

ER^{−} | 10 | 2.8 ± 2.2 |

Progesterone receptor (PR) status | ||

PR^{+} | 84 | 9.7 ± 8.6 |

PR^{−} | 27 | 5 ± 3.9 |

HER2 status | ||

HER2^{+} | 11 | 8 ± 5.8 |

HER2^{−} | 94 | 8.5 ± 7.6 |

HER2 equivocal | 6 | 11.7 ± 15.5 |

Lymph node status | ||

Positive | 54 | 9.1 ± 7.1 |

Negative | 49 | 9.2 ± 9.15 |

Unknown | 8 | 3.2 ± 1.7 |

Total | 111 carcinomas (of 251 samples) |

Diagnosis/type . | Number . | |${\overline {| {\rm {G^*}} |}}$| (kPa) . |
---|---|---|

Histological type | ||

Invasive ductal carcinoma (IDC) | 75 | 8.8 ± 8.5 |

Invasive lobular carcinoma (ILC) | 15 | 10.5 ± 5.3 |

Invasive ductal lobular carcinoma (IDLC) | 7 | 5.9 ± 3.7 |

Invasive mucinous carcinoma (IMC) | 3 | 4.4 ± 1.4 |

Invasive papillary carcinoma (IPC) | 8 | 2.4 ± 1.3 |

Invasive micropapillary carcinoma (IMPC) | 3 | 19.7 ± 12.6 |

Histological grade | ||

Grade 1 | 10 | 13.9 ± 12.2 |

Grade 2 | 46 | 10.2 ± 8.7 |

Grade 3 | 55 | 6.3 ± 5.3 |

Estrogen receptor (ER) status | ||

ER^{+} | 101 | 9.2 ± 8.1 |

ER^{−} | 10 | 2.8 ± 2.2 |

Progesterone receptor (PR) status | ||

PR^{+} | 84 | 9.7 ± 8.6 |

PR^{−} | 27 | 5 ± 3.9 |

HER2 status | ||

HER2^{+} | 11 | 8 ± 5.8 |

HER2^{−} | 94 | 8.5 ± 7.6 |

HER2 equivocal | 6 | 11.7 ± 15.5 |

Lymph node status | ||

Positive | 54 | 9.1 ± 7.1 |

Negative | 49 | 9.2 ± 9.15 |

Unknown | 8 | 3.2 ± 1.7 |

Total | 111 carcinomas (of 251 samples) |

### Mechanical rheometry

A conventional strain-controlled rheometer (ARG2, TA instruments Inc.) was used to evaluate |G*(ω)| for a subset of *N* = 28 breast tissue specimens (Supplementary Data S5). The sample was sandwiched between the rheometer plates (8 mm diameter) and a frequency sweep was performed over ω = 0.628 to 250 rad/second while maintaining the strain below 0.5%, to obtain G*(ω).

### Statistical analysis

Linear regression analysis investigated the correlations between SHEAR measurements, collagen content, and rheometer measurements. Two-way ANOVA was used to evaluate associations between prognostic criteria and SHEAR measurements (Prism, GraphPad Inc.). Pairwise comparisons were performed using Tukey method (for multiple comparisons), unpaired (between independent groups) or paired *t* tests (between dependent groups). Multivariate analysis and a generalized estimating equation model were used to identify prognostic criteria that significantly and independently correlated with |$\overline {\ | {{G^*}} |}$| and |$\overline {\ | \nabla |{G^*}||} $| (SPSS, IBM Corp.). In all cases, *P* < 0.05 was considered statistically significant.

## Results

### SHEAR |G*(x,y)| maps correspond with breast histopathology

Figure 2A shows the |G*(x,y)*|* map of a normal breast specimen. Within adipose regions (H&E image; Fig. 2B), low |G*(x,y)*|* values approximately 0.05 to 3 kPa were observed (white arrows), whereas in fibrous tissue |G*(x,y)| increased to 5 to 20 kPa, consistent with the bright pink hues in H&E image and the collagen fiber network in PSR slide (Fig. 2C, magenta arrows). In glandular tissue, including ducts and terminal duct lobular units (TDLU), |G*(x,y)| was approximately 3 to 5 kPa (yellow and red arrows), with higher moduli observed near ducts than around TDLU, consistent with the expected stiffness variations in the vicinity of glandular tissues (22). Linear regression analysis indicated a significant positive correlation between the spatially averaged log( |$\overline {| {{{\rm{G}}^{\rm{*}}}} |} $| ) obtained by SHEAR and collagen content, measured as the percentage area occupied by orange–yellow hues of circularly polarized light microscopy images of PSR-stained sections (*N* = 251, *r* = 0.69, *P* < 0.0001; Fig. 2D). The *y*-axis represented the log( |$\overline {| {{{\rm{G}}^{{*}}}} |}$| ), because while collagen density saturates at 100%, fiber alignment and cross-linking could substantially increase |$\overline {| {{{\rm{G}}^{\rm{*}}}} |}$| even for small changes in total collagen.

### Frequency-dependent SHEAR measurements correlate with standard rheometry

In Fig. 2E–G, the spatially averaged modulus, |$\overline {| {{{\rm{G}}^{\rm{*}}}( {\rm{\omega }} )} |} $| measured by SHEAR is compared with that measured by a mechanical rheometer in 28 breast tumors at three distinct frequencies. At low frequencies, |$\overline {| {{{\rm{G}}^{\rm{*}}}( {\rm{\omega }} )} |} $| measured by SHEAR exhibits a strong correlation with corresponding rheometer measurements (at ω = π rad/second: *R* = 0.78, *P* < 0.0001; and ω = 2π rad/second: *R* = 0.76, *P* < 0.0001) over a large range of moduli values (0.7–17 kPa). We further observe a good correlation with rheometer measurements, which is maintained up to ω = 10π rad/second (5Hz, *R* = 0.64, *P* = 0.0002). At ω = 20π rad/second, the correlation drops to *R* = 0.41 (*P* = 0.04) likely due to inertial effects of the rheometer at higher frequencies (Supplementary Data S5).

### SHEAR maps are influenced by histopathological subtypes of breast cancer

#### |${\overline {| {{G^{\rm{*}}}} |}} $| is excessively elevated in aggressive invasive micropapillary carcinoma compared with other subtypes

Figure 3A, E, I, and M present the |G*(x,y)| maps of typical invasive mucinous carcinoma (IMC), invasive ductal carcinoma (IDC), invasive papillary carcinoma (IPC), and invasive micropapillary carcinoma (IMPC) tumors, the corresponding |∇|G*(x,y)|| (Fig. 3B, F, J, and N), H&E images (Fig. 3C, G, K, and O), and PSR sections (Fig. 3D, H, L, and P). In IMC (Fig. 3A), regions that exhibit |G*(x,y)|∼0.2 to 3 kPa coincide with the mucinous compartments in the H&E image. These compliant areas contrast the regions of increased |G*(x,y)| ∼30 kPa, corresponding to collagen-rich stroma in the PSR image. In the |G*(x,y)| map of Fig. 3E, which represents an IDC tumor with a ductal carcinoma *in situ* (DCIS) component, the tumor epithelium exhibits low values of |G*(x,y)| ∼0.1 to 5 kPa in regions where tumor cells are confined within ducts (DCIS, yellow arrows) and when tumor cells breach out as invasive carcinoma (red arrows). Moreover, the necrotic cores of DCIS exhibits a lower |G*(x,y)| ∼3 kPa. Unlike epithelial regions, |G*(x,y)| exceeds 5 kPa within the stroma: around DCIS, |G*(x,y)| is ∼10 kPa, whereas for desmoplastic stroma of the IDC component, |G*(x,y)| ∼10 to 20 kPa. The moderately aggressive IPC subtype in Fig. 3I exhibits distinct features with very low |G*(x,y)| ∼0.1 to 1 kPa corresponding to the highly vascularized epithelium, comprising low viscosity blood, and a central stromal stem of slightly increased |G*(x,y)| ∼5 kPa. Figure 3M displays the |G*(x,y)| map of IMPC, which is associated with a particularly high rate of metastasis and poor prognosis. The reduced |G*(x,y)| ∼3 kPa areas correspond to highly malignant cells in the H&E image with thin spikes of increased |G*(x,y)| that extend from the stroma and closely correspond to collagen architecture in the PSR section. These bundled and stiffened collagen fibers, aligned perpendicular to IF, are believed to support the preferential migration and invasion of tumor cells (12). A highly contrasting |G*(x,y)| ∼30 kPa in Fig. 3M coincides with the accumulation, thickening, and occasional perpendicular realignments of collagen fibers in the PSR image of Fig. 3P. Consistent with |G*(x,y)| maps of Fig. 3, ANOVA of 251 human breast specimens from 148 patients in our cohort reveals that |$\overline {| {{{\rm{G}}^{\rm{*}}}} |} $| is significantly different among benign and tumor specimens of distinct histological subtype (*P* < 0.0001; Fig. 3Q). Tukey multiple comparisons indicate that adipose tissue exhibits a significantly lower |$\overline {| {{{\rm{G}}^{\rm{*}}}} |} $| than benign fibrous, IDC, invasive lobular carcinoma (ILC), and IMPC. In addition, |$\overline {| {{{\rm{G}}^{\rm{*}}}} |} $| in both fibrous tissues and IPC is significantly lower than IMPC, given the prominently desmoplastic stroma of the latter. Multivariate analysis indicates that histologic subtype is significantly associated with |$\overline {| {{{\rm{G}}^{\rm{*}}}} |} \ $| (*P* < 0.004) with aggressive IMPC tumors exhibiting highly elevated |$\overline {| {{{\rm{G}}^{\rm{*}}}} |} \ $| (Supplementary Tables S1 and S2).

#### The |∇|G*(x,y)|| is significantly increased in aggressive subtypes

The IF interface between the epithelium and stroma are delineated by an abrupt increase in |G*(x,y)| and intensified |∇|G*(x,y)|| for all tumors in Fig. 3. Among all lesions, the |∇|G*(x,y)|| of IMPC in Fig. 3N exhibits the highest value at the IF (white arrow). ANOVA performed in 38 of 111 tumor specimens that present a clear IF based on the H&E, suggests that the |$\overline {| \nabla |{G^*}||} $| at IF is significantly different among subtypes, being highest in IMPC and lowest in IPC (Fig. 3R). Subsequently, multivariate analysis confirms that histologic subtype is the most influential prognostic criteria independently associated with |$\overline {| \nabla |{\rm G^*}||} $| (*P* < 0.0001; Supplementary Tables S3 and S4).

### |G*(x,y)| maps of breast carcinoma are influenced by tumor grade

#### |${\overline {|{\bf G^{\rm{*}}}|}} $| is reduced in high-grade tumors

Figure 4A, E, and I show |G*(x,y)| maps of tumors with increasing grades, the corresponding |∇|G*(x,y)|| (Fig. 4B, F, and J), H&E images (Fig. 4C, G, and K), and PSR sections (Fig. 4D, H, and L). The |G*(x,y)| within the tumor epithelium reduces with grade from approximately 6.5 kPa in grade 1 (Fig. 4A), to approximately 4.8 kPa in grade 2 (Fig. 4E), and approximately 2 kPa in grade 3 (Fig. 4I) tumors. In contrast, the |G*(x,y)| of the tumor stroma, in grade 1 and 2 tumors with perceptible stromal components, is considerably elevated to 30 kPa as seen in the left half of Fig. 4A and bottom half of Fig. 4E, consistent with densely packed collagen fibers seen in PSR sections of Fig. 4D and H. In the grade 2 tumor of Fig. 4E, sporadic bands of increased |G*(x,y)| of up to 7 kPa in the compliant regions correspond to sparse intratumoral ECM fibers in the epithelium, as evidenced by the H&E image (Fig. 4G). Unlike Fig. 4A and E, the |G*(x,y)| map of the grade 3 lesion in Fig. 4I is mostly compliant and the H&E image (Fig. 4K) reveals a highly cellular tumor. Moreover, the extremely reduced |G*(x,y)| < 2 kPa in Fig. 4I correspond to mitotic cells. These contrast the neighboring areas of |G*(x,y)| ∼5 kPa (red arrows) that coincide with necrosed cells (yellow arrows) in Fig. 4K. From Fig. 4, it appears that the amount of peripheral stroma, and in turn the overall |$\overline {| {{\rm G^*}} |} $| , is reduced in high-grade lesions. ANOVA consequently confirms that |$\overline {| {{\rm G^*}} |} $| reduces with tumor grade, and grade 3 lesions exhibit the lowest |$\overline {|{\rm G^*}|} $| values (Fig. 4M; *P* < 0.05). Multivariate analysis subsequently indicates that tumor grade is significantly associated with |$\overline {| {{\rm G^*}} |} $| (*P* < 0.002; Supplementary Tables S1 and S2).

#### High-grade tumors exhibit increased |$\overline {| \nabla |{{\rm{G}}^{\rm{*}}}||} $| at the IF

Strong mechanical contrast between the tumor epithelium and stroma leads to an increased |$\overline {| \nabla |{{\rm{G}}^{\rm{*}}}||} $| at the IF, as evidenced by Fig. 4B and F. In contrast, grade 3 tumor depicted in Fig. 4I is more necrotic than fibrotic and lacks a perceptible IF. Nevertheless, a subset of grade 3 lesions in our cohort, for example, Fig. 3M, exhibit considerable peritumoral fibrosis and an intensified |$\overline {| \nabla |{{\rm{G}}^{\rm{*}}}||} $| at the IF. ANOVA results show that |$\overline {| \nabla |{{\rm{G}}^{\rm{*}}}||} $| is significantly different among tumors of different grades, with grade 3 lesions exhibiting larger gradients than grade 2 tumors at the IF (Fig. 4N). However, when considered collectively with other prognostic criteria in a multivariate analysis, tumor grade is not independently associated with |$\overline {| \nabla |{{\rm{G}}^{\rm{*}}}||} $| .

#### |$\overline {|{\bf G^{\rm{*}}}|} {\rm{\ }}$| is significantly lower in the tumor epithelium versus stroma

To investigate heterogenous moduli distributions, |G*(x,y)| maps were separated into epithelium and stroma components (Fig. 5A, E, and G), through multiplication with spatial masks (Fig. 5D; Supplementary Data S4). This analysis was performed in 58 (of 111) tumors in which both components could be clearly seen in H&E images (Fig. 5B, F, and H). Figure 5A, E, and G show distinct variations in |G*(x,y)| distributed within epithelium and stroma of a grade 1, IDC tumor: significantly lower moduli (∼10 kPa) are observed within the epithelium versus the stroma (∼20–30 kPa).

The corresponding |∇|G*(x,y)|| is considerably elevated at the epithelium–stroma interface (Fig. 5C). The pooled data (Fig. 5I and J) of |$\overline {| {{G^*}} |} $| indicate that both epithelial and stromal moduli significantly reduce with the tumor grade (*P* = 0.0001). Paired *t* tests further shows that within each grade, stroma is significantly stiffer than epithelium. Moreover, grade 1 tumors entail both a stiffer stroma and epithelium than grade 3 lesions, suggesting that malignant, high-grade cells are progressively more compliant (33).

Intratumoral heterogeneities of |G*(x,y)| maps were also estimated by calculating the information Entropy (Supplementary Data S4). Entropy was measured in 106 (of 111) tumors, for which the |G*(x,y)| maps were free of mosaicking artifacts. We observe that Entropy trends toward lower values with increasing tumor grade (Fig. 5K), supporting our above findings that the |G*(x,y)| in high-grade cellular masses is uniformly low.

### Mechanical characteristics of tumor microenvironment are driven by receptor expressions

#### Molecular expressions are linked to elevated |${\overline {| {{G^{\rm{*}}}} |}}$|

Figure 6A, E, and I show |G*(x,y)| maps of grade 3 tumors with different molecular expressions, the corresponding |∇|G*(x,y)|| (Fig. 6B, F, and J), H&E images (Fig. 6C, G, and K), and PSR sections (Fig. 6D, H, and L). Epithelial areas exhibit lower |G*(x,y)| (∼0.1–5 kPa), where tumor cells mingle with adipose tissue, degenerated ECM fibers, and immune infiltrates. In contrast, the stroma exhibits distinctively higher |G*(x,y)| values. In the ER/PR^{+}, HER2^{−} tumor, |G*(x,y)| increases up to 20 kPa (Fig. 6A), corresponding with bright red collagen bundles (Fig. 6D) within the fibrotic stroma. In ER/PR^{+}, HER2^{+} tumor, |G*(x,y)| features islands of moderately increased |G*(x,y)| up to 10 kPa, harboring spots of elevated |G*(x,y)|∼30 kPa (Fig. 6E) that correspond with thick collagen fibers and microcalcifications in H&E. Finally, the ER/PR^{−}, HER2^{−} (triple-negative breast cancer, TNBC) lacks a significant fibrotic component (Fig. 6K and L) and presents moderate values of |G*(x,y)| around necrosed areas. From these |G*(x,y)| maps, it appears that TNBC tumors exhibits a significantly reduced |$\overline {| {{\rm G^*}} |} $| compared with ER/PR^{+}, likely due to higher cell proliferation and more necrosis than fibrosis in this aggressive subtype (34). Subsequently, unpaired *t* tests indicates that |$\overline {| {{\rm G^*}} |} $| is significantly higher in ER^{+} and PR^{+} lesions (Fig. 6M and N). Multivariate analysis, however, shows that only the PR status significantly affects |$\overline {| {{\rm G^*}} |} $| (*P* = 0.043; Supplementary Tables S1 and S2). No significant association is identified in |$\overline {| {{\rm G^*}} |} $| with respect to HER2 status (Fig. 6O).

#### |$\overline {\nabla | {{\bf G^{\rm{*}}}} |} {\rm{\ }}$| is increased at the IF of tumors expressing ER, PR, and HER2 receptors

The |∇|G*(x,y)|| is elevated at the IF of the ER^{+}/PR^{+}/HER2^{−} lesion (Fig. 6B), whereas regions of elevated |∇|G*(x,y)|| are observed in the vicinity of microcalcifications for the ER^{+}/PR^{+}/HER2^{+} tumor (Fig. 6E). In the TNBC tumor (Fig. 6I), |∇|G*(x,y)|| is mildly elevated around necrotic regions. Subsequently, unpaired *t* test shows that |$\overline {| \nabla |{{\rm{G}}^{\rm{*}}}||} $| is significantly higher in ER^{+} and PR^{+} tumors (Fig. 6P and Q). For these tumors, a higher |$\overline {| {{\rm G^*}} |} $| is noticed as well, likely dominated by the stromal components, which in turn give rise to higher |$\overline {| \nabla |{\rm G^*}||} $| at the IF. The |$\overline {| \nabla |{{\rm{G}}^{\rm{*}}}||} $| also trends higher in HER2^{+} lesions, although this is not significant (Fig. 6R). Multivariate analysis indicates that only ER status (*P* = 0.023) independently affects |$\overline {| \nabla |{{\rm{G}}^{\rm{*}}}||} $| (Supplementary Tables S3 and S4).

### Micromechanical characteristics of tumors with lymph node involvement

SHEAR maps displayed in Figs. 3A and M, 4E, 5A, and 6A are from tumors with lymph node involvement (LN^{+}). Although majority of LN^{+} specimens exhibit a trend toward increased |G*(x,y)| in tumor periphery, multivariate analysis does not find a significant association in |$\overline {| {{{\rm{G}}^{\rm{*}}}} |} $| with respect to LN status (Supplementary Fig. S4A). In contrast, |$\overline {| \nabla |{\rm G^*}||} \ $| is significantly higher in LN^{+} specimens (*P* < 0.0001; Supplementary Fig. S4SB). Multivariate analysis show that LN status (*P* = 0.047) is independently associated with |$\overline {| \nabla |{{\rm{G}}^{\rm{*}}}||} \ $| (Supplementary Tables S3 and S4). These findings agree with previous hypotheses on directional invasion of cancer cells, suggesting that the metastatic migration is derived by the |$\ \overline {| \nabla |{{\rm{G}}^{\rm{*}}}||} $| , rather than overall tumor |$\overline {| {{{\rm{G}}^{\rm{*}}}} |} $| (35, 36).

## Discussion

Breast carcinoma are frequently classified on the basis of histological subtype, grade, molecular receptors’ expression, and lymph node involvement. These prognostic criteria likely both influence and are influenced by the viscoelastic modulus, |G*(x,y,ω)| of the tumor microenvironment (1–3, 7, 12, 21, 22, 37). In clinical settings, ultrasound and magnetic resonance elastography (MRE) have been shown to detect the increased elasticity of tumor tissue (38–42). Nevertheless, due to low resolution, these bulk-scale measurements remain insufficient for assessing nuanced associations of tumor histopathology with stiffness. Recently, atomic force microscopy (AFM) enabled high-resolution, albeit laborious, mapping of subsurface elastic properties across small FoVs (∼10s of μm), over hours, in frozen sections and demonstrated the link between micromechanical heterogeneities and invasive subtypes (36, 43). Moreover, optical coherence elastography enabled tumor margin delineation based on microscale deformations of breast tissue over FoVs of multiple cm^{2} (44, 45). However, both techniques contact the tissue to apply a static load and require *a priori* knowledge of stress to evaluate the strain response, and thus only indirectly assess the static elastic properties, but fail to recapitulate the dynamic viscoelastic behavior of tumors (36, 43–45).

Here, we have demonstrated for the first time the capability of SHEAR to estimate and map the frequency-dependent viscoelastic modulus, |G*(ω)|, with microns-scale resolution, over centimeters of FoVs in excised human breast specimens, in a noncontact manner. SHEAR captures multiply-scattered, multi-speckle patterns; therefore, speckle intensity fluctuations are sensitive to even 10′s of nm displacements of individual scattering particles. This in turn permits probing |G*(ω)| over a wide range of highly compliant to very stiff breast tissue specimens (Fig. 1; Supplementary Fig. S1). Using homogenous hydrogel materials, we have previously demonstrated that |G*(ω)| values measured from speckle intensity fluctuations were closely correlated with both rheometer and AFM measurements and spanned a large moduli range of 0.1 to 40 kPa (29). We have extended these prior studies to demonstrate the strong correlation between |G*(ω)| measured by SHEAR and mechanical rheometry in the breast tissue specimens at multiple frequencies (Fig. 2E–G). The larger spread of data points, compared with our prior studies in homogenous hydrogels, is likely due to challenges associated with conducting rheometry in tissue, differences between the length scales of rheometer and SHEAR measurements, and tissue heterogeneity (Supplementary Data S5). The range of moduli measured in our study are comparable with those previously measured by MRE in typical benign and malignant breast tissue (42), however, SHEAR measurements were significantly larger than previous AFM results. This was in part due to the smaller length scale and FoV, and superficial probing depth of AFM that likely weighted the measured elastic moduli toward lower values in these prior studies.

To estimate |G*(ω)| via GSER, the sphere-equivalent particle radius, *a*, was estimated by comparing g_{2}(t) at perpendicular and parallel polarization states, g_{2⊥}(t) and g_{2||}(t), respectively (Supplementary Data S2). A limitation of our study is that a single value of *a* = 100 nm was used for all 251 specimens based on the above methods. However, scattering particles’ radius, *a*, could vary both among different tissue types and within different regions of tissue. For example, in a representative set of 23 breast tissues, we estimated that *a* varied by up to ±35 nm, suggesting a maximum variability of 35% in the estimated |G*(ω)| values, which likely contributed in part to the spread of |G*(ω)| in the scatter plots of Fig. 2E–G. In future, this variability may be significantly reduced by calculating spatially varying ratios of g_{2⊥}(t) and g_{2||}(t) and estimating distinct particle radii while scanning the sample.

In addition to scattering particle size, biological tissue present a variety of length scales, with the ECM mesh structure varying over a large range of pore sizes from tens of nm to a few microns (46). Therefore, depending on the relative size of scattering particles, *a*, and pore size, ξ, the MSDs of scattering particles probe distinct mechanical behaviors of the tissue constituents (47). More specifically, in weaker, unlinked ECM regions, where *a*<ζ, the motion of scattering particles is largely diffusive and the MSD reflects the local viscoelasticity of the cellular, necrotic, or lipid components depleted of ECM linkages. On the other hand, in dense ECM network, where ξ<*a*, the scattering particles are effectively coupled to the matrix and the MSD probes the viscoelastic properties of the cross-linked ECM network. This in fact is precisely the source of high mechanical contrast in SHEAR measurements. In other words, in the context of heterogeneous tissue, pixel-by-pixel analysis of multi-speckle frames in SHEAR enables probing local scattering particles whether they are restricted in small pores of dense and desmoplastic ECM, or diffuse between the looser, un–cross-linked meshwork. Thus, replacing the MSD and *a* in the GSER yields the local viscoelastic properties, at scales larger than the size of scattering particles, yet small and inaccessible to macroscopic rheometry (47, 48). This could likely further explain the spread in the scatter plots of Fig. 2E–G.

The frequency range of |G*(ω)| is mainly governed by the acquisition duration and frame rate (Supplementary Data S1). The Basler acA2000-340 km camera permitted a frequency limit of 6,000 rad/second; however, in our study we limited the frame rate to 250 fps with the goal of comparing SHEAR results with mechanical rheometry, which is limited by tool inertia at higher frequencies.

While this study is limited to the investigation of the shear modulus, |G*(ω)|, additional information may be gained by extracting the relative contributions of G'(ω) and G"(ω). In an ongoing pilot study, we observed that in highly desmoplastic tumors, with restricted particle displacements, G' dominates the G*, whereas for highly cellular, vascular, or mucinous lesions, with diffusive particle displacement, G" is likely the prominent contributor. While representative curves are presented in this article (Fig. 1; Supplementary Fig. S1) to demonstrate the opportunity to estimate G'(ω) and G"(ω), a comprehensive investigation is yet necessary to fully investigate the complex associations between distinct elastic and viscous moduli with tumor prognostic criteria. Similarly, for ease of interpretation with histopathology, shear moduli maps at a single frequency ω = 2π rad/second are presented. Yet, SHEAR may offer a wealth of information to investigate whether the frequency-dependent evolution of |G*(x,y,ω)| presents additional, unique micromechanical markers for clinical prognostication.

Our results further confirmed the significant differences in average |$\overline {| {{{\rm{G}}^{\rm{*}}}} |} \ $| among tumors of different prognostic criteria. The overall |$\overline {| {{{\rm{G}}^{\rm{*}}}} |} $| was higher in aggressive histologic subtypes, and in tumors enriched in hormonal and HER2 receptors, but reduced with histologic grade, and remained indifferent with respect to lymph node status. The increased |$\overline {| {{{\rm{G}}^{\rm{*}}}} |} $| in ER^{+} specimen agreed with previous accounts on enriched ER expression in mammographically dense patients (49). In addition, the reduced |$\overline {| {{{\rm{G}}^{\rm{*}}}} |} $| in TNBC lesions was corroborated by prior ultrasound elastography studies and attributed to the elevated necrosis and diminished desmoplastic reaction in these lesions (50).

The distinct differences in |G*(x,y)| within the tumor epithelium and stroma components provided a unique perspective on the intratumoral heterogeneity. Intriguingly, certain micromechanical features were observed in aggressive tumors. For instance, the tumor epithelium in higher grade lesions was highly compliant, a finding that is corroborated by previous studies that suggest acute deformability of extremely malignant tumor cells (8, 33). Furthermore, our finding that tumor stroma is significantly stiffer versus epithelium across all grades is supported by the well-established concept of the desmoplastic reaction (2). In addition, the reduction of entropy with grade likely suggests that in high-grade masses, compliant, highly proliferative cells yield a mechanically uniform landscape. Increased |$\overline {| {\nabla | {{\rm G^*}} |} |} $| at the IF of ER^{+}, PR^{+}, HER2^{+}, and LN^{+} lesions was observed in our study. The role of stiffness gradients in directing local invasion has been demonstrated *in vitro* by using engineered gels of tunable elasticities and tracking the migratory behavior of seeded cells (6). Through characterizing the |$\overline {| {\nabla | {{\rm G^*}} |} |} $| in clinical specimens, SHEAR may facilitate the translation of these fundamental mechanobiology concepts toward clinical diagnosis. By evaluating and identifying the micromechanical features shared between aggressive subgroups, SHEAR likely provides a unified understanding of the invasive pathophysiology beyond traditional phenotyping criteria (51).

## Authors' Disclosures

Z. Hajjarian reports grants from American Society for Laser Medicine and Surgery and Eleanor and Miles Shore Faculty Development Awards Program during the conduct of the study, as well as a patent for compensation for causes of temporal fluctuations of backscattered speckle patterns in laser speckle rheology of biological fluids issued, a patent for system and methods estimation of mechanical properties and size of light-scattering particles in materials issued, a patent for optical thromboelastography systems and methods issued, a patent for laser speckle micro-rheology in characterization of biomechanical properties of tissues issued, and a patent for compensating for variations of optical properties in laser speckle microrheology of breast lesions pending. D.M. Tshikudi reports a patent for optical thromboelastography systems and methods issued. S.K. Nadkarni reports grants from NIH during the conduct of the study; grants and nonfinancial support from Coalesenz Inc. outside the submitted work; a patent for optical thromboelastography system and method for evaluation of blood coagulation metrics pending, issued, and licensed to Coalesenz Inc., a patent for evaluation of viscoelastic properties of tissue based on laser speckle fluctuations issued, a patent for holographic speckle microrheometer (HSM) for measuring micromechanic micromechanical properties issued, a patent for compensation for multiple scattering in laser speckle rheology of biological fluid issued, a patent for device for comprehensive assessment of blood coagulation in real time issued and licensed to Coalesenz Inc., a patent for system and methods for estimation of mechanical properties and size of light-scattering properties in materials issued, a patent for optical thromboelastography systems and methods issued and licensed to Coalesenz Inc., a patent for laser speckle micro-rheology in characterization of biomechanical properties of tissues issued, a patent for apparatus, devices and methods for obtaining omnidirectional viewing by a catheter pending and issued, a patent for non-contact optical system, computer-accessible medium and method for measuring at least one mechanical property of tissue using coherent speckle issued and licensed, and a patent for compensating for variations of optical properties in laser speckle microrheology of breast lesions pending; and holds equity in and serves on the board of directors of Coalesenz Inc., which is developing a point-of-care system to measure blood coagulation parameters. S.K. Nadkarni's interests were reviewed and are managed by Massachusetts General Hospital and Partners Healthcare in accordance with their conflict-of-interest policies. No disclosures were reported by the other authors.

## Authors' Contributions

**Z. Hajjarian:** Conceptualization, data curation, software, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing. **E.F. Brachtel:** Formal analysis, writing–review and editing. **D.M. Tshikudi:** Data curation, methodology, writing–review and editing. **S.K. Nadkarni:** Conceptualization, resources, data curation, software, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing.

## Acknowledgments

The authors thank Arriane Garret for technical assistance. They sincerely thank Mohmmad Miri, Maxwell Hogue, Marina Kem, Rea Ruiz, Carolyn Mc Donagh, Haley Martin, and Tim Casson from MGH tissue repository for providing the human breast tissue specimens. The authors also thank Hang Lee of the Harvard Catalyst Program for providing statistical consultation on data analysis.

This work was supported in part by the NIH [NHLBI 5R01HL119867, principal investigator (PI): S.K. Nadkarni], an American Society for Laser Medicine and Surgery research grant (PI: Z. Hajjarian), and an Eleanore and Miles Shore Foundation Fellowship (PI: Z. Hajjarian).

The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.

**Note:** Supplementary data for this article are available at Cancer Research Online (http://cancerres.aacrjournals.org/).