Multiparametric Classification Links Tumor Microenvironments with Tumor Cell Phenotype

Tumor microenvironment features are established as predictors of tumor cell behavior and fate.


Introduction
Statistics from 2013 show that more than 90% of breast cancer deaths are a consequence of metastatic disease and that median survival for metastatic breast cancer is 3 years [1]. Accordingly, a better understanding of the physiological conditions that help initiate or facilitate metastatic recurrence would be valuable for identifying and improving outcomes for patients with metastatic breast cancer. The process of metastasis [2] includes epithelialmesenchymal transition, tumor cell locomotion through the interstitial matrix, dissemination from the primary tumor via vasculature (intravasation and extravasation), and seeding of secondary organs [3]. The chemical signals that initiate and direct tumor cell motility are typically growth factors secreted either by host cells in the primary tumor microenvironment (e.g., macrophages, endothelial cells, fibroblasts, and pericytes) or by tumor cells themselves [4], and may be present either in soluble form or bound to the extracellular matrix (ECM). Recently, increasing evidence has shown that mechanical signals in tumor microenvironments may also modify tumor cell behavior. Tumor cell locomotion has been shown to correlate with proximity to macrophages [5] and blood vessels [6] while increased crosslinking [7], stiffness [8,9], and alignment [10] of the ECM correlate with increased metastasis.
In summary, many studies have investigated individual tumor cell-stromal cell or tumor cell-ECM interactions, revealing a number of possible sources of chemo-and hapto-tactic signals for tumor cells, as well as the high level of complexity in the tumor microenvironment network [4]. Moreover, a number of tumor microenvironment studies [11] have shown positive or negative effects of stromal cell types or individual mechanical parameters on tumor cell motility within the primary tumor and subsequent metastasis. However, an integrated view of the tumor microenvironment as a complex system, including how all relevant biological players in combination affect tumor cell motility and metastatic outcome, is still missing.
In ECM that is soft, porous, and with low levels of cross-linking, single tumor cells can locomote quickly and in a metalloprotease (MMP)-independent manner [12] by utilizing forces generated by adhesion and actin polymerization [13,14]. Locomotory protrusions (pseudopodia) form at the front of the cell, followed by locomotion and translocation of the entire cell. Stiff, noncompliant ECM contains high concentrations of cross-linked collagen and/or basement membrane proteins. To enable movement, tumor cells remodel their immediate ECM, mainly via the degradative action of MMPs [12,15,16]. Studies in cell culture models have shown that, when grown on ECM in the presence of growth factors, metastatic cancer cells form invadopodia [17]. These are dynamic membrane protrusions that are enriched in actin, actin regulatory proteins, cortactin, tyrosine kinases, Tks5, and proteases [18], and are ECM-degrading structures. In vitro, invadopodia form in 10%-20% of tumor cells plated on ECM [19], and their average size in the light microscope is approximately 1-3 mm in diameter and 2-12 mm in length, depending on ECM density and dimensionality (2-D versus 3-D) [19][20][21].
The invadopodium compartment is simultaneously the site of actin polymerization, cell contact to ECM, and matrix proteolysis; thus, the invadopodial lifetime is on the scale of hours, significantly longer than any other type of membrane protrusions in tumor cells [19,22]. Although a large part of the invadopodial signaling network overlaps with that of other membrane structures such as lamellipodia, filopodia, and focal adhesions, the defining feature of invadopodia is their high proteolytic activity, which degrades and remodels ECM. This suggests that invadopodia may be a key mechanism for MMP-dependent tumor cell invasion [23].
Recent studies have gone beyond 2-D and into the 3-D realm of studying invadopodia [20] and podosomes [24]. Moreover, there are recent attempts to study these structures in situ, in their natural environment, including podosomes in mouse aorta [25], invadopodia in dermis [26], as well as protrusions similar to invadopodia in organogenesis models of Caenorhabditis elegans [27] and C. intestinalis [28]. Further, we have recently shown a direct link between the ability of tumor cells to assemble invadopodia and degrade matrix using MMPs in vitro and in situ and the ability to intravasate and metastasize in vivo [29,30]. Collective data from these and similar studies demonstrate that morphology and protein distributions of invadopodia and podosomes depend on the matrix dimensionality, architecture, and complexity as well as the cell line/type. Neither tumor cells nor the structures they assemble appear the same in 2-D or 3-D [20] culture, or within the fixed tissue [29]. Hence, to identify podosomes and/or invadopodia in different microenvironments and under new experimental conditions, a combination of defining features such as small size, structural components (e.g., Tks5, cortactin) and functionality (MMP-dependent ECM-degradation) may be used, classifying structures as podosomes when referring to structures in myelomonocytic, endothelial, and vascular smooth muscle cells [17] and invadopodia when referring to tumor cells.
Although tumor cell locomotion in primary tumors has been observed in a number of studies [4,5,[31][32][33][34], the exact location and timing of invadopodium formation is not presently known. Combined results from numerous studies using preparations of collagen I (with or without cross-linking), collagen IV, laminin, fibrin, chorioallantoic, and peritoneal membranes [12,35] converge on the conclusion that invadopodia have the ability to degrade and remodel different ECM components present in interstitial tissue and basement membranes of metastatic tumor models [29]. However, where and when invadopodia form in the primary tumor in vivo has yet to be determined. Here, we have characterized and quantified two motility phenotypes occurring in primary breast tumors in vivo: fast-locomotion and slow-locomotion associated with invadopodia. We show that invadopodia (and slow-locomotion) occur in regions of primary tumor that are spatially distinct from regions where fastlocomotion is the dominant phenotype. We further combined intravital multiphoton microscopy, image analysis, and multiparametric data classification, to analyze microenvironmental conditions under which tumor cells move via either phenotype. We found that the prediction of primary tumor locations where tumor cells locomote fast or form invadopodia and locomote slowly cannot be achieved by any individual tumor microenvironment parameters but only by their combination.

Results
By recording 4-D stacks in orthotopic xenograft tumors formed from a metastatic, breast cancer cell line MDA-MB-231-Dendra2, we confirmed the presence of a tumor cell behavior previously described in detail, where cells quickly locomote through the tissue either individually or in streams ( Figure 1A-1C; Movie S1a)

Author Summary
A large proportion of cancer deaths are due to metastasis-the spread of cancer from the primary tumor to other parts of the body. Movement of cells may require the formation of protrusions called invadopodia, which degrade extracellular matrix. Although some studies have reported on locomotion in primary tumors, the presence of invadopodia was not tested. Here, we show that single cells from mouse mammary carcinoma can move using a fast-or slow-locomotion mode depending on different levels of cues present in the tumor microenvironment. Using multiphoton microscopy in vivo combined with a machine-learning algorithm we show how manipulation of microenvironmental conditions can induce predictable changes in the number of locomoting cells or switch between the two locomotion modes. We also demonstrate that only the slower moving cells are associated with invadopodia in vivo, and that only tumor cells from regions rich in invadopodia degrade the surrounding extracellular matrix and disseminate. Specific targeting of invadopodia results in inhibition of lung metastasis. This work proposes a systems biology view of how tumor microenvironments regulate tumor progression and presents insight into the heterogeneity of the treatment response. The ability to define and predict conditions under which tumor cells disseminate offers potential therapeutic benefits in regulating tumor progression. [4,[31][32][33][34]. Subtraction of the initial frame (09) from that taken at 30 min (309) in a single z-section, detects the parts of the cell that moved over 30 min and the area the cell covered (D30 cell tracks) ( Figure 1B, blue). However, in some fields of view, tumor cells exhibited a different phenotype. Fast-locomotion was not detected ( Figure 1D and 1E; Movie S1b) but was replaced by slowlocomotion. Imaging at high magnification revealed the presence of small protrusions ( Figure 1E, white arrowheads and red cell tracks; Movie S2), often adjacent to collagen fibers (Movie S2a and S2b), macrophages and blood vessels (Movie S2c and S2d). Longer time-lapse imaging showed that tumor cells with small protrusions exhibited a slower locomotion phenotype Figure 1F and 1G, red) in comparison to previously reported fast-locomotion ( Figure 1F and 1G, blue). Cells migrating with velocities of 22-250 mm/h were classified as fast-locomoting, while those with small protrusions and velocities 2-15 mm/h were classified as slowlocomoting ( Figure 1F).
Systematic analysis of fields of view in ten animals has shown that the two motility phenotypes were usually mutually exclusive, i.e., tumor cells within a specific field of view exhibited only one of the observed phenotypes at a time in 184/187 fields ( Figure 2A). As xenograft tumors consist of genetically similar tumor cells, we hypothesized that the presence of two motility phenotypes is controlled via the tumor microenvironment and may be distinguished by different levels of one or more microenvironmental factors, whose measurement can be then used as a predictor of tumor cell behavior.
To investigate which microenvironmental conditions are amenable for either motility phenotype, we monitored tumor cell phenotype relative to the number of microenvironment parameters previously reported to influence tumor cell locomotion ( Figure 2B) [5][6][7][8][9][10]. These included the density of collagen fibers [7][8][9][10], tumor cells, and macrophages [5] as well as the number and size of blood vessels [6] present in the field of view. We tested each of the parameters individually for their correlation with either of the motility phenotypes, attempting to find values of a single parameter or combination of them that coincide with one of the two phenotypes. However, none of the commonly used statistical methods based on linear separation (model selection, partial correlation, or causal modeling) (Figure S1A-S1C) were able to provide a model capable of distinguishing between microenvironment parameters and the presence of different motility phenotypes. A possible reason for the lack of linear predictors is the existence of complex interactions among parameters, such as those present in the tumor microenvironment [4,11]. To test this, we employed support vector machine (SVM), an algorithm suitable for segmentation of non-linearly separable data, such as the data in Figure S1A. Intuitively, the SVM algorithm ''converts'' a nonlinearly separable dataset into linearly separable set by increasing the dimensionality of the hyperspace (initially, N parameters create hyperspace of N dimensions). For illustrative purposes, see a onedimensional example in Figure S11A.
In order to link between microenvironmental parameters and tumor cell locomotion phenotypes, we used microenvironmental parameters as an ''input'' for the SVM classifier, with fastlocomotion (blue dots) or slow-locomotion (red dots) as the classifier's ''output'' ( Figure S1D). Classification was done on the basis of three, four, or five parameters with increasing accuracy, reaching maximum accuracy of 92% when all five parameters were used ( Figure S1D). Figure 2C shows a 3-D projection of the hyperspace with maximum accuracy classification. The high accuracy of the SVM classification (and cross-validation) means that it is possible to predict whether tumor cells will locomote fast or slow (output) on the basis of a given set of microenvironmental parameters (input).
By utilizing only the microenvironmental parameters within a field of view as input information, we get a prediction of the motility phenotype present in that same field of view, which is output information. Misclassification (green dots) probability of ,8% (and ,5% when training includes the entire dataset) is associated mainly with parameter values on the classification border between red space and blue space ( Figure 2C, red-and blue-shaded areas). None of these individual microenvironment parameters alone offers a sufficient predictor, as their contribution to motility phenotype depends on the context created by all other microenvironment parameters. Only a nonlinear transformation of all parameters could distinguish between microenvironments associated with either motility phenotype ( Figure S1). In conclusion, tumor motility phenotypes can be distinguished only by a nonlinear, multiparametric analysis such as SVM, as they are a result of balance among multiple signals within complex tumor microenvironments.
To better understand the link between the tumor microenvironment and tumor cell motility phenotypes, we analyzed each of the motility phenotypes at single cell level ( Figure 3). We characterized the relationship between individual tumor cells of either phenotype and landmarks within the tumor microenvironment ( Figure 3A). While directionality of fast-locomotion seemed to be controlled by the orientation of adjacent collagen fibers ( Figure 3A, blue pie charts), the small protrusions in slowlocomoting cells were commonly perpendicular to adjacent blood vessels and collagen fibers that surround blood vessels ( Figure 3A,  pie charts). To further analyze the relationship between blood vessels and small protrusions, we compared the distance of tumor cells with either motility phenotype to the nearest blood vessel and the size of that blood vessel ( Figure 3B). Our results show that fastlocomotion may occur far from blood vessels (21/132 of fastlocomoting cells were found.100 mm from the nearest blood vessel) and is independent of the blood vessel size. However, the distances from the nearest vessel at which small protrusions formed were smaller (2/95 of slow-locomoting cells were found.100 mm from the nearest blood vessel) and positively correlated with vessel diameter ( Figure 3B, red triangle), suggesting that small protrusions are most likely initiated by a chemical or mechanical gradient associated with the blood vessel. We also measured the amount of ECM components not detected by second harmonic generation ( Figure S2), demonstrating that deposits of basement membrane proteins laminin ( Figure S2A) and collagen IV ( Figure S2B) increase with blood vessel size, suggesting that small protrusions are associated with larger blood vessels with this matrix composition.
We next looked at the number of motile cells within fields of view where either fast-or slow-locomotion was detected. Both fast (blue bars) and slow (red bars) locomotion each occurred in approximately 15% of tumor cells in the field of view ( Figure 3C) and only slow-locomoting cells exhibited small protrusions ( Figure 3D). Finally, small protrusions associated with slowlocomotion were found to be inhibited by tail-vein injection of MMP inhibitor GM6001, suggesting MMP-dependence, while the fast-locomotion remained unaffected throughout the 3 h time lapse ( Figure 3E).
Our phenomenological data, including morphology, relationship to ECM and blood vessels, as well as MMP-dependence of small protrusions led us to the hypothesis that the observed small protrusions are invadopodia in vivo. To test this hypothesis, we used previously established structural and functional markers of invadopodia. Molecularly, invadopodia have been shown to be enriched in active proteases and structural components including cortactin and Tks5 [35][36][37][38][39]. While cortactin is present both in invadopodia and at the leading edge of migrating cells in 2-D culture [40], in 3-D conditions it is enriched primarily at the tip of invadopodial protrusions at the cell front [20,29]. By using MDA-MB-231-cortactin-GFP cells [20], we were able to directly compare cortactin-enriched compartments in 3-D collagen I ( Figure S3A), in cryosections ( Figure S3B) [29] and in vivo ( Figure 4). Similarly to 3-D and cryosections, the small protrusions observed in vivo (left and middle panels) showed a peak of cortactin fluorescence at the protrusion tip ( Figure 4A, yellow lines in upper panels and associated line-scans in lower panels; Movie S4b). In contrast, fast-locomoting cells showed a homogeneous distribution of cortactin throughout the cell ( Figure 4A, right panels; Movie S4a). These results are consistent with the identification of the small protrusions as invadopodia in vivo.
We next constructed cell lines where the Tks5 structural protein of invadopodia was knockdown (KD) with shRNA, which specifically targets invadopodia in MDA-MB-231 cells and tumors ( Figure S4). Tks5 knockdown was previously shown to inhibit invadopodium maturation and degradation [39,41] in 2-D culture conditions [39] as well as to greatly reduce lung metastases in mouse models [41]. Our data confirmed that Tks5 KD1 and Tks5 KD2 cells ( Figure S4A) did not assemble invadopodia in vitro, under culture conditions ( Figure S4B and S4C), as shown using Tks5/cortactin antibodies and fluorescent gelatin degradation. Similar results were confirmed in primary tumors, where knockdown efficiency was maintained ( Figure S4D), and Tks5positive punctae which colocalized with collagen I degradation ( Figure S4E) were present only in TKs5 CTRL tumors but not in knockdowns ( Figure S4F). In addition, Tks5 KD1 and Tks5 KD2 tumors in vivo did not exhibit small protrusions, while the fast locomotion behavior was only slightly affected ( Figure S4G), supporting our hypothesis that small protrusions are indeed invadopodia and they were selectively targeted by Tks5 knockdown.
Finally, to directly test if the small protrusions function as invadopodia in vivo, we measured the ECM degradation activity associated with the small protrusions. The antibody against degraded collagen we used for immunohistofluorescence ( Figure  S3) is unsuitable for in vivo use due to the inefficient delivery and labeling. Instead, we used the MMP-activated substrate MMPSense 680 (Perkin Elmer) for intravital imaging [42]. To validate this reporter, we compared ECM degradation as measured by MMPSense 680 solution (cyan) and a more commonly used substrate, DQ-collagen I gel (red) [14,43] in 3- Quantitation of ECM degradation area with or without MMP inhibitor GM6001 ( Figure S5C) showed similar trends with both reporters. In vivo, MMPSense 680 (cyan) injected into MDA-MB-231-cortactin-GFP tumors (green) colocalized with cortactinenriched protrusions after 3 h ( Figure 4B, left panels; yellow lines correspond to line-scans; Movie S5a). In contrast, fast-locomoting cells did not colocalize with MMPSense 680 ( Figure 4B, right panels; Movie S5b). Finally, we monitored the accumulation of MMPSense 680 signal in microenvironments where small protrusions were either present or absent, in microenvironments treated with GM6001, as well as in Tks5 knockdown tumors. This experiment was done using photoconversion followed by MMPSense 680 injection and image collection from 0-96 h ( Figure S6). We quantified MMPSense 680 signal ( Figure 4C) at 24 h, when the signal in the tumor plateaus. In fields of view where small protrusions were present, we measured that 6% of image area contained MMPSense 680 fluorescence above threshold, reflecting the locations of MMP-based degradation and this was approximately 3-fold higher than in other areas, GM6001-treated areas or knockdown tumors. In addition, MMPSense 680 fluorescence was positively correlated with the number of small protrusions ( Figure 4D), and this correlation was eliminated by GM6001 treatment. As all our results are consistent with the hypothesis that the small protrusions are in fact invadopodia in vivo, we will henceforth refer to them as invadopodia.
The results of SVM classification show that the two motility phenotypes, which we now recognized as fast-locomotion and invadopodia-associated slow-locomotion, exist in distinct yet neighboring conditions. Such a result suggests that a shift along one of the axes in the 3-D plot in Figure 2C may induce a change in the number of tumor cells exhibiting either motility phenotype, i.e., induce a phenotype switch. The changes in tumor cell behavior are context-dependent, which means that the extent of the influence of one parameter over the phenotype, or its capability to cause a phenotype switch, depends on the context of other parameters ( Figure S7). We hypothesized that a motility phenotype switch in primary tumors may occur, for example, when collagen fiber density changes over time. An increase in collagen fiber density may be initiated because of fibroblast deposition of collagen, while decreases may be a result of the invadopodium-mediated degradation of fibers. Such logic is strengthened by previous in vitro reports showing that both the speed of MMP-dependent 3-D migration [44] and the number of invadopodia in 2-D assays are controlled by the rigidity and crosslinking level in basement membrane extracts, collagen, and synthetic matrices [45,46]. We tested the effect of ECM rigidity/ cross-linking by modulating ECM cross-linking levels and measuring the number of invadopodia, which are associated with slowlocomotion phenotype ( Figure 5). In the control set of animals, we imaged the same fields of view (using photoconversion to match fields over time) at 0, 24, and 48 h, demonstrating that invadopodia are present over the entire period under control conditions ( Figure 5A). Collagen imaging confirmed that over a 48 h period, collagen I fibers remained stable with minor changes ( Figures 5E, purple bars, and S8A). A different set of animals was treated with L-ribose, which was shown to increase cross-linking and hence stiffness in collagen-based gels [7]. A considerable increase in the number of invadopodia ( Figure 5B) was detected, accompanied by an increase in collagen I density (Figures 5E, green bars, and S8B). The third set of animals was treated with the lysil oxidase (LOX) inhibitor b-aminopropionitrile (BAPN), which reduces cross-linking and loosens the ECM [7]. Treatment with BAPN reduced and finally eliminated invadopodia ( Figure 5C). Figure 5D shows that invadopodium number follows the trends in collagen fiber density. Moreover, as a result of fibrillar collagen I density decrease (Figures 5E, red bars, and S8C), in several fields of view, tumor cells switched from invadopodium-associated slowlocomotion at 0 h to fast-locomotion at 48 h in the same field of view ( Figure 5C, inset).
Areas with fast-locomotion also showed a dramatic reduction of fast-locomoting cells following the L-ribose treatment ( Figure 5F, green bars) accompanied by the appearance of invadopodia. Interestingly, BAPN treatments also induced a slight decrease in the number of fast-locomoting cells over 48 h ( Figure 5F, red bars). One of the possible explanations may be the lack of balance between adhesion and traction forces at low ECM stiffness levels, which was previously shown to affect tumor cell migration in 3-D environments [44].
The observation of phenotypic switching is consistent with the prediction of SVM classification ( Figure 2C), which although none of the microenvironment parameters individually are sufficient predictors of tumor cell motility phenotype, a change in one or more of them in the context of the others can modulate and even switch the tumor cell behavior. Therefore, a decrease in collagen fiber density, depending on the context of other parameters ( Figure S7) may transform the microenvironment that favors the invadopodial phenotype ( Figure 2C, red sphere) into the microenvironment that favors fast-locomotion ( Figure 2C, blue sphere), in turn switching the motility phenotype; similarly, an increase in collagen fibers may induce a switch from fast-to slow-locomotion is worth noting that a similar phenotypic switch, from MMPdependent to MMP-independent motility, has been previously demonstrated in vitro by inhibiting ROCK or MMPs in melanoma cells [47,48] and fibrosarcoma cells [49]. However, previous to our results reported here, it was unclear if such a phenotypic switch can occur in vivo and what collagen density would regulate this phenotypic switch.
Finally, we hypothesized that the two motility phenotypes may lead to different tumor cell fates ( Figure 6). To test this, we photoconverted subpopulations of tumor cells within microenvironments where either slow-( Figure 6A) or fast-locomotion ( Figure 6B) was present. A 1756175 mm 2 square was photoconverted within each of the four to 12 neighboring fields of view at 0 h (see Figure S6), enabling us to follow the tumor cell fate at subsequent timepoints (Materials and Methods and [6,29] for experimental details). In the fields with invadopodia, there was a significant negative trend in the number of photoconverted (red) cells at 24-48 h ( Figure 6C, red bars), previously shown to be due to dissemination from the primary tumor [6,29,32,50]. In contrast, other fields that contained fast-locomoting cells had cell numbers that were mostly unchanged and sometimes increased in numbers of photoconverted cells, likely a result of the cell division ( Figure 6C, blue bars) [6]. In addition, the number of invadopodia at 0 h shows very strong negative correlation (p = 5.4610 26 ) with the number of photoconverted cells at 24 h ( Figure 6D), suggesting that the presence of invadopodia is directly linked to tumor cell disappearance from the field of view and, possibly, tumor cell dissemination [3]. Next, GM6001 irreversibly eliminated invadopodia ( Figure 6E) and inhibited the disappearance of red cells ( Figure 6C, grey bars), supporting the hypothesis that tumor cells with invadopodia disseminate from the primary tumor [29].
To confirm that the disappearance of red cells in the areas rich with invadopodia (but not in areas with fast-locomoting cells) is due to the intravasation, we tested alternative processes that may contribute to changes of red tumor cell numbers, such as apoptosis and failure to proliferate (Figures S9 and S10). By combining photoconversion in vivo with additional ex vivo labeling ( Figure  S9), it is possible to asses rates of proliferation ( Figure S9A and S9C) and apoptosis ( Figure S9B and S9C) in fields of view with either fast-or slow-locomoting cells. No significant differences were observed in the proliferation rates of fields of view associated with slow-and fast-locomotion and there was practically no apoptosis observed. To confirm such an observation, an injectable marker SR-FLIVO was used for direct apoptosis measurements in  real time ( Figure S10), combined with cisplatin treatment to induce apoptosis in the positive control group. This method also reported the absence of apoptosis in areas where fast-or slowlocomotion occurs, strengthening the link between the red cells disappearance and intravasation. Following disseminated red cells to the secondary organs, we monitored the number of tumor cells that had formed lung metastases at 5 days post-photoconversion ( Figure 6F and 6G). Some lung metastases contained red tumor cells (photoconverted cells that arrived at the lung at 0-5 days from the photoconversion site in the primary tumor), orange and yellow tumor cells (photoconverted cells where red/green Den-dra2 ratio is reduced due to the cell division and synthesis of new green Dendra2), and green tumor cells (whose origin cannot be traced) ( Figure 6F). When tumors were chemically treated with BAPN (which selectively inhibits invadopodium formation by modifying the external microenvironment) or GM6001 (which inhibits ECM degradation by invadopodia) starting at photoconversion time (0 h, see Materials and Methods), we observed a significant reduction in the number of metastases containing red tumor cells ( Figure 6G, panel i), pointing to onset of inhibition of dissemination at 0 h. Moreover, lung metastases were practically eliminated in Tks5-knockdown tumors ( Figure 6G, panel ii) where invadopodia were removed, while the fast-locomotion phenotype remained largely unaffected ( Figure S4G). These data together support the hypotheses that invadopodia and ECM degradation by tumor cells are essential for dissemination and metastasis, and that metastasis does not occur without active ECM degradation by tumor cells.

Discussion
In this study, we monitor two different phenotypes of tumor cell motility: fast-locomotion and invadopodium-associated slow-locomotion, which appear in spatially separate microenvironments in primary, orthotopic xenograft tumors. Such a separation of phenotypes among different fields was striking and suggests that the size of the imaging field is well below the average size of each phenotypic microenvironment. As both phenotypes can be seen in the same animal (imaging window covers ,50 mm 2 ) but they did not commonly occur in the same (,0.5 mm 2 ) or adjacent fields, we estimate the phenotypic microenvironment size to be 1-50 mm 2 .
Although both motility phenotypes are dependent on the presence of vasculature, ECM, and macrophages, they do not show obvious dependency on any of these parameters, either individually or when (linearly) combined. However, by using a multiparametric, non-linear SVM classification, we were able to predictably classify microenvironmental parameters that favor each motility phenotype. The non-linear nature of this statistical model exposes the complexity of tumor microenvironment, a theme often addressed in reviews, but commonly avoided in experimental science, where the focus is on individual microenvironmental parameters and regard others as stable and homogeneous ''background''. Here we show that more than one microenvironmental component is involved in phenotype determination in primary tumors and also, that small regions with genetically identical tumor cells may exhibit different phenotypes. We have analyzed five microenvironment parameters (number of macrophages, density of collagen fibers and tumor cells, number and size of blood vessels), which all contribute to the determination of tumor cell motility phenotype to a smaller or larger extent. The extent of contribution of each parameter depends on the context of the other four parameters ( Figure S7) and can change through time. Hence, heterogeneity, complexity, and dynamical changes in the tumor microenvironment dominate tumor cell phenotype and should be considered in the interpretation of future intravital imaging studies and in the development of new diagnostics and treatments. Stepping away from the focus on single microenvironment parameter analysis towards the multiparametric analysis may not only contribute to generating a more comprehensive view of how tumor microenvironments determine cell behavior, but also help with the search for new therapeutic targets, which is currently conducted by analysis of individual molecules or specific protein interactions.
SVM classification suggests the existence of a fine balance among chemical and mechanical signals produced by tumor cells, macrophages, blood vessels, and ECM, which determines the motility phenotype of tumor cells. Ongoing changes in one or more microenvironment parameters (collagen fibrillogenesis, angiogenesis, macrophage recruitment, etc.) can lead to a switch between fast-locomotion and invadopodium-associated slow-locomotion. It is likely that both motility phenotypes contribute to metastasis. Fast-locomotion is more efficient in translocation of tumor cells, which may transport them to regions where they stop due to ECM composition and stiffness, such as in regions adjacent to some blood vessels as shown here, and then shift to the invadopodium-associated phenotype that results in tumor cell dissemination. The presence of misclassifications ( Figure 2C, green dots) on the borders between two phenotypes suggests that there may be microenvironmental conditions where the switch between motility phenotypes commonly takes place. In further support of this, we were able to modulate invadopodium frequency by changing one of the parameters, collagen fiber density ( Figure  S6). A decrease in collagen fiber density resulted in elimination of invadopodia and, in some areas, a switch to fast-locomotion. In theory, changing other parameters such as blood vessel size and macrophage number would achieve a similar result ( Figure S6). To understand the mechanism of this phenotypic switching, further mathematical modeling and direct in vitro measurements are currently underway. We speculate that both motility phenotypes are likely to be exhibited by the same tumor cells and that there are no genetic mutations unique to the two motility phenotypes. In support of this, we see both phenotypes in a xenograft mouse model, based on genetically identical tumor cells. Interestingly, the fraction of cells exhibiting either of the phenotypes in a particular field of view is the same, 15%. Similar fractions of tumor cells were previously reported as capable of assembling invadopodia in situ, tumor sections, or in vitro, when plated on thin gelatin [29]. This leads to the suggestion that both fast-and slow-locomotion are phenotypes that arise when similar internal conditions are present (motility cycle is active) but the cell is exposed to different external conditions. We propose that within the same region, under the same external conditions, the other 85% cells are not capable for motility at the moment, and the reason for that may be that they are mitotic or bound by energetic and signaling constrains. Broadening this hypothesis into other motility phenotypes, similar fractions of cells were also reported in individually or collectively migrating cells in mammary carcinoma xenografts, where phenotype was dependent on intracellular transforming growth factor beta (TGFb) expression [51].
We demonstrated that the slow-locomotion phenotype, which was not previously characterized in vivo, is associated with small tumor cell protrusions identified as invadopodia using morphological, molecular, and functional assays. This motility phenotype is linked to ECM degradation, tumor cell dissemination from the primary tumor and subsequent lung metastasis, and inhibited by MMP-inhibitor. It is worth noting that the treatment with GM6001 was previously shown in 2-D culture to minimally affect existing invadopodial core structures [52], but their growth and extension were not tested. In fact, invadopodia in 3-D grow in cycles of extension-retraction [20] and degradation of ECM by MMPs followed by enlargement or elongation of invadopodia was suggested to be a part of a positive feedback mechanism [37]. We propose that by abolishing this cycle, GM6001 inhibits elongation, rendering protrusions too small to be detected in vivo. In summary, our results support the hypothesis of a positive feedback loop between products of MMP proteolysis and continued invadopodial extension [37]. Based on the evidence collected in vitro in 2-D and 3-D [12,15,53], we hypothesize that membranetethered MT1-MMP1 is likely enriched in the actively degrading invadopodia compared to fast-locomoting cells. However, none of the available antibodies were able to give us a satisfactory signal. While experimentally, MMPs remained only broadly tested using GM6001, we speculate that the key role of MMPs includes immobilization of MT1-MMP on invadopodia, leading to local activation of all MMPs present in the external microenvironment [54].

Ethics Statement
All procedures involving animals were conducted in accordance with NIH regulations, and approved by the Albert Einstein College of Medicine IACUC, by 20101010 and 20130909 protocol numbers.

Cell Lines and Animal Models
A metastatic human breast cancer line MDA-MB-231 was cultured and maintained in DMEM media supplemented with 10% fetal bovine serum (FBS) and 50 U penicillin/50 mg streptomycin per ml. The Dendra2-MDA-MB-231 cell line [50] was created by electroporation of the Dendra2 cloning vector C1 with the geneticin selection marker [6]. We performed fluorescence-activated cell sorting (FACS) after 14-20 days of selection; after removing the highest-expressing top 5%, we kept the top 20% of the highest fluorescing cells and maintained them under selection using 500 mg/ml geneticin (Invitrogen). No changes in cell morphology, viability, or proliferation were observed in the labeled MDA-MB-231 cells compared to the parental cell line. Cell lines labeled with cortactin-GFP [20] and cortactin-TagRFP [19] have been previously described. Tks5 knockdown cell lines Tks5 KD1 and Tks5 KD2 and knockdown control cell line Tks5 CTRL were all created by transduction of Dendra2-MDA-MB-231 cells with lentiviral particles (five particles/cell), which contained shRNA in pLKO.1 vector, targeting Tks5 (knockdowns) or non-targeting shRNA (CTRL) (Sigma Aldrich MISSION library based on Broad Institute consortium). Tumors were formed by injecting 10 6 cells in 150 ml of 20% collagen I in PBS into the mammary fat pads of 5-7 week old SCID mice. Imaging experiments were done 9-12 weeks after injection.

Multiphoton Microscopy
Multiphoton imaging was done on a custom-built system based on an inverted Olympus IX71 microscope and two Ti-Sapphire lasers, one of them equipped with an Optical Parametric Oscillator (OPO) extension [55]. The first laser was tuned to 880 nm for excitation of GFP-like fluorophores and Texas Red, while the second laser was tuned for photoconverted Dendra2 at 1035 nm and MMPSense 680 at 1250 nm [55]. For experiments involving simultaneous Dendra2 and MMPSense 680 imaging, MMPSense680 was excited using the 880 nm laser. The microscope system is equipped with four simultaneously acquiring detectors, which allowed simultaneous imaging of collagen fibers via second harmonic generation (blue) and fluorescence from GFP-like fluorophores (green), Texas Red or photoconverted Dendra 2 (red), and MMPSense (far red).
In the photoconversion-enabled fate mapping experiments, a commercially available microscope system with multiphoton and confocal modes (Olympus FV1000MPE with ULTRA 256, 1.05 NA water immersion objective) was used for initial 4-D imaging over 30 min (in multiphoton mode, with laser locked at 880 nm) and subsequent photoconversion of Dendra2 (in confocal mode, using 405-, 488-, and 546-nm laser lines).

Continuous Imaging Sessions Using Skin Flap
Multiphoton imaging of mammary tumors was done as described previously [56]. MDA-MB-231 cells labeled with Dendra2-or cortactin-GFP were injected into the mammary fat pads of 5-7 week old SCID mice. After 10-12 weeks, we injected Texas Red 70 kDa dextran (250 nmol in 100 ml PBS per injection) [5] for macrophage labeling and performed skin flap surgery 2 h later on anesthetized animals. Exposed mammary tumors were positioned on top of a coverslip on an inverted microscope and imaged continuously using a custom-built two-laser multiphoton microscope for up to 3-5 h. In some experiments, additional tailvein injections were done for blood vessel labeling using Texas Red 70 kDa dextran (Invitrogen; 250 nmol in 100 ml PBS per injection), MMPSense 680 (Perkin Elmer; 2 nmol in 150 ml PBS per injection), SR-FLIVO (Immunochemistry Technologies; 1:10 dilution, 2 h prior to imaging), or the pan-MMP inhibitor GM6001 (Milipore; 1 mmol in 100 ml injection). A stock solution of GM6001, 500 mM in DMSO, was diluted in sterile PBS before tail vein injection. Post-surgical injection assures that compounds are only present in tumor blood vessels that are intact, connected to tumor vasculature and flowing at the time of injection.

Photoconversion-Enabled Tumor Cell Fate Mapping through the Mammary Imaging Window
As previously described [6,29,50,[57][58][59], the mammary imaging window was combined with photoconvertible Dendra2 as a cytoplasmic marker. Briefly, Dendra2-MDA-MB-231 cells were injected into the mammary fat pads of 5-7 week old SCID mice. After 9-11 weeks, a mammary imaging window was implanted on top of the growing tumor (5-7 mm in diameter). Following a 3-day recovery, mice were anesthetized and placed in the imaging chamber on an Olympus FV1000-MPE hybrid multiphotonconfocal system. In each area, ten z-sections were acquired (total stack size is 512 mm6512 mm6100 mm) over 30 min. One to three areas were imaged per animal. Subsequently, photoconversion was done in 1756175 mm areas, separated by 150 mm along both x and y axes ( Figure S6). Photoconversion areas were scanned for ten to 20 cycles with a 405-nm laser in confocal mode. Further imaging was done using the two-laser multiphoton system, using 880 nm for imaging of green species of Dendra2 and 1,035 nm for red species of Dendra2. 3-D stacks were collected at photoconverted sites at 0 h, 24 h, and 48 h post-photoconversion. The number of red (photoconverted) Dendra2 cells was counted from the RGB overlay on individual z-sections, using manual Cell Counter plugin in ImageJ. Reported numbers are sums of counts at five z-sections separated by 25 mm, for method details see [29]. Values were corrected for cell division by dividing by control values from areas where no motility was detected [6]. The dispersion of photoconverted cells is monitored as the surface embedding all red cells in the 0-100 mm maximum image projection [6].
In MMPSense 680-based degradation measurements (Figure 4B-4D), in SR-FLIVO apoptosis measurements ( Figure S9), and in microenvironment modulation experiments ( Figure 5), photoconversion was used to locate the same microenvironments at 0-48 h. Intravenous injections of 150 ml MMPSense 680 were done at 0 h, intraperitoneal injections of 150 ml PBS containing either 1 g/kg/day L-ribose or 20 mg/kg/day BAPN were done at 0 h and 24 h, while intraperitoneal injections of 2 mg/kg/day cisplatin were done daily for 5 days.

Photoconversion-Enabled Tumor Cell Fate Mapping in Pulmonary Metastases
As previously described [55], transdermal photoconversion of Dendra2 can be used for fate mapping in pulmonary metastases. Briefly, 12 weeks after injection of Dendra2-MDA-MB-231 cells, pairs of SCID mice were matched by tumor size from the same cage. In each pair, one mouse was injected intraperitoneally with 150 ml of sterile PBS (control), while the other was injected with 1 mg GM6001 in 150 ml PBS (GM6001-treated [31]). Six hours later, mice were anesthetized by isoflurane and hair was shaved off the top of the mammary tumor. The tumors were photoconverted transdermally, using a 405-nm LED photodiode array [55]. Treatments with PBS or GM6001 were repeated daily, and mice were sacrificed at 5 days. Lungs were taken out and metastatic colonies (groups containing 1+ tumor cells) were counted in 50 fields of view through the ocular of the microscope, using a 256 objective. Colonies containing 1+ red tumor cells were classified as red, although most also contained green cells owing to the cell division over 5 days. Representative images ( Figure 6F and 6G) were taken using an Olympus multiphoton setup, using the same imaging parameters as for intravital imaging. Counted colonies were analyzed by contingency table and compared for their average values ( Figure 6H), resulting in significant differences between control and GM6001-treated red metastases (x 2 = 20.04, p = 0.048).
Labeling of tumor cell proliferation and apoptosis postphotoconversion was done by cryosectioning of the primary tumor regions under the mammary imaging window (top 150 mm, [60]). Neighboring sections were labeled by Ki67 (DAKO) and cleaved caspase 3 (Milipore), using Alexa633 as secondary antibodies.

Blood Vessel Coverage by ECM Components
Blood vessel coverage by ECM components, collagen I and IV, laminin and fibronectin was measured using ImageJ ( Figure S2). Briefly, images of blood vessels and ECM proteins were thresholded and tested for colocalization using Manders' M2 coefficient [61]. Blood vessel edges were traced using the green channel, with a 5-mm-thick band generated to encapsulate the area of ECM proteins. The amount of ECM components deposited around blood vessels was measured and normalized to the surface area.

MMPSense 680 Validation
MMPSense 680 (Perkin Elmer) is a polymer with quenched VT680 chromophore and is weakly fluorescent prior to exposure to MMPs. Proteolytic cleavage of the peptide that links the copolymer probe can be executed by a wide range of MMPs and produces a signal in the near-infrared range (with emission peak at 680 nm) [62]. As MMPSense 680 is soluble, it can be added either to culture media or injected directly into the animal vasculature.
To validate MMPSense 680 as a reporter of MMP-driven proteolytic activity at single cell resolution ( Figure S5), a comparison was made to an established marker of matrix degradation, FITC-DQ-collagen I (Invitrogen) in 3-D culture based on rat-tail collagen I (BD Biosciences). A 2-4-mm-thick layer of DQ-collagen I (10 mg/ml) and collagen I (2 mg/ml) mixture (1:100, DQ-unlabeled collagen) in cold, phenylalanine-free DMEM/10% FBS serum was established on the bottom of a Mattek dish. After 20 min at 37uC, 5610 4 MDA-MB-231cortactin-TagRFP cells were seeded on this layer and covered with additional collagen I. After 60 min at 37uC, we added 50 nM MMPSense 680 and cultured the cells for 6-24 h. Image stacks were collected at 0, 6, and 24 h using a Leica SP5 confocal microscope, utilizing 488-, 543-, and 633-nm laser lines and setting the prism for detection at 500-525, 555-610, and 670-740 nm. Image processing and quantification of the amount of degraded ECM, measured as percent of total image area after threshold, was done in ImageJ. Briefly, 3-D stacks of degraded DQ collagen or MMPSense 680 were combined into a maximum intensity projection and processed using a 1-pixel median filter. Thresholding of images taken at 24 h ( Figure S5A) was done using average values of images taken at 0 h.

Intravital Systems Microscopy
The workflow of our intravital systems microscopy approach is outlined in Figure 2A and 2B. Each 4-D image stack was collected in three or four channels. After a brief quality screening of each zsection, two adjacent z-sections were combined into a maximum intensity projection (10 mm thick) and the recorded stack was separated into four neighboring fields of view (2566256 mm) for easier processing. All images were passed through smoothing filters and thresholded to remove background fluorescence. Individual microenvironment time lapses were passed through the StackReg/ TurboReg plugin [63], available at http://bigwww.epfl.ch/ thevenaz/stackreg/. These Java classes efficiently remove minor drifts in the xy plane resulting from breathing and tissue settling.
Time-lapse movies of individual z-sections in each field of view (309-609 duration) were first visually scored for morphological determinants of tumor cell fast locomotion, as previously described [4][5][6]29,[31][32][33][34]. Briefly, the fast locomoting cell translocates by the extension of the cell front, movement of the center and the contraction of the rear, with an average speed of 1 mm [34]. To account for the number of fast locomoting cells the frame taken at 09 was subtracted from the frame taken at 309, both in the green channel, yielding an image of all pixels translocated during this time interval. Particle analysis and Region of Interest (ROI) Manager in ImageJ were then used to count number of locomoting cells and overlay them with the original 09 frame to show initial tumor cell position and the trace left by the tumor cell migration over 309 ( Figure 1B, green for tumor cells and blue for the tumor cell tracks). Some tumor cells moved out of the analyzed z-section during the 309 period and were traced in the neighboring z-sections. In such cases, cell tracks were visualized via maximum projection, accounting for the thickness of the z-slice (5 mm).
Next, time-lapse movies were viewed at 2-46 zoom, exposing small protrusions characterized by (a) finger-like morphology, (b) width of 1-3 mm, and (c) dynamics. Frame subtraction and particle analysis were used to visualize the dynamics of small protrusions over 309 (tip movement and change in size due to the cycles of extension and retraction) as well as the number of small protrusions per cell ( Figures 1F, red line, and 2D). In order to detect migration of the tumor cell body (front, center, and rear) in cells with small protrusions, movement of time-lapse movies were extended to 5 h. The following five microenvironment parameters were extracted from the multichannel 3-D stacks: (a) density of collagen fibers, defined as percent area above threshold in the blue channel, (b) density of tumor cells, defined as percent area in the green channel, (c) number of macrophages, defined as macrophages labeled with 70 kDa dextran [5], (d) number of flowing blood vessels in the field of view, and (e) diameter of the largest flowing blood vessel visible in the field of view. For (d) and (e), we made the assumption that all blood vessels flowing at the time of post-surgical injection by Texas Red were labeled. As macrophages and blood vessels were collected in the red channel, prior to measurements they were separated on the basis of size and morphology, using the Analyze Particles plugin in ImageJ.
Microenvironment parameters and the motility phenotype (i.e., fast-or slow-locomotion) form a multiparametric matrix that we used for SVM classification [64]. Different types of kernels and other parameters were tested systematically in R software. Classification was finally done using Gaussian radial basis kernel, using the ''tune.svm'' procedure in the R-package ''e1071,'' starting from the default cand cost-parameters and iteratively optimizing them to increase accuracy, which we tested using 10fold cross-validation. The dataset was divided into randomly grouped ten equal-size sets, training on nine while testing on the left-out set. Starting from SVM based on three-parameter combinations, we were able to achieve 78%-84% accuracy ( Figure S1D). By adding all five parameters measured, we ultimately achieved 92% accuracy. Without cross-validation (i.e., using all data) classification accuracy is 95%. This means that we were able to predict motility phenotype on the basis of measurements of all five microenvironment parameters in 92/ 100 fields of view and that removing any of the parameters reduces the predictive power of SVM. However, for purposes of presentation in 3-D space, we have used a 3-D projection that contains parameters with the highest influence: density of collagen fibers, number of macrophages, and diameter of blood vessels ( Figure 2C). Misclassifications ( Figure 2C, 5/100, green dots) are present in boundary areas. Note that blue-and red-shaded areas in Figure 2C are provided only as an illustration.

Data Processing and Analysis
Collected image stacks were processed by a combination of the ImageJ [65] and Imaris software packages. ImageJ was supplemented with custom-written plugins based on our work and work of others for image browsing, color correction, and migration trajectories. We used Imaris 7.4 for 3-D reconstructions. Data was directly exported to Microsoft Excel for statistical analysis and plotting; the R package was used for trend analyses, variance trend analysis, SVM classifications, correlations, linear regression, and calculation of p-values, while Matlab was used for SVM representation. All numerical values used in figure graphs are included in Data S1 file.

Variance Trend Analysis
To test the influence of the blood vessel distance and diameter on tumor cell motility phenotypes ( Figure 3B), we measured the shortest distance of the tumor cell front to the nearest blood vessel. We hypothesized that in case the blood vessel is the source of a chemical or mechanical gradient, then the increase in blood vessel diameter will allow motility at larger distances, thus increasing the spread of distances (variance) of motile tumor cells from the blood vessel. To test this hypothesis, we first calculated the residuals that result from regressing the observed distance of the cell front on the vessel diameter. In turn, the log square residuals were regressed on the vessel diameter. As hypothesized, an increase in blood vessel diameter was associated with a significant increase in residuals in cells with small protrusions (R 2 = 0.3, p = 7610 28 ) indicating that in the fields of view where only microvessels are present, small protrusions form only next to the blood vessel, while in the presence of large-diameter blood vessels, small protrusions can be both close and far from the blood vessel. Fast-locomoting cells exhibit no such behavior (R 2 = 0.07, p = 10 23 ). Note the small R 2 value. Red-shaded area in Figure 3B illustrates variance trend as blood vessel diameter increases. Figure S1 Distinguishing between tumor microenvironment conditions where fast-or slow-locomotion occurs can not be done with linear models but with non-linear SVM modeling. (A) Model selection using any of 1, 2, 3, 4, or 5 microenvironment parameters does not provide separability of microenvironments in which fast locomotion (blue) or slowlocomotion (red) occur. As an illustration, all 2-parameter scatter plots are shown. There are no lines that can be drawn to separate blue and red points without significant misclassification. N MF, Number of macrophages; N BV, Number of blood vessels; D BV, diameter of blood vessels; % ECM, area covered by collagen fibers; % TC, area covered by tumor cells. (B, C) Partial correlation or causal modeling do not yield statistically significant results. (D) Non-linear classification using SVM methodology resulted in 92% accuracy, with only a few misclassified points. Accuracy is based on 10-fold cross-validation. Presented here are two SVM classifications based on three out of the five microenvironmental parameters. On the left is a classification based on D BV, N MF, and % ECM, yielding 84% accuracy; on the right, classification is based on N BV,% ECM, and % TC, yielding 78% accuracy. Misclassified data (green points) are tumor cells whose observed class of protrusion is opposite to that predicted by the SVM classification algorithm. Numerical values used for analysis are the same as in Figure 2C. A shift from level 2X to level 1X in the parameter X will cause a phenotypic switch from slow-locomotion (red) to fast locomotion (blue) motility when parameter Y is at level 1Y, but not when parameter Y is at level 2Y. In the experiment shown in Figure 5, parameter X is collagen fiber density, which we manipulated by changing ECM cross-linking using ribose or BAPN. Parameter Y may be any of the other four parameters (number of macrophages or blood vessels, blood vessel diameter, or tumor cell crowding). Data S1 Numerical data used to generate plots in Figures 1-6 and S1, S2, S3, S4, S5, S6, S7, S8, S9, S10, S11.  Figure 1E). Examples of small protrusions directed towards collagen fibers (panels a, b) or blood vessels (panels c, d). Some small protrusions in direct contact with blood vessel are protruding into the blood vessel (panels c, d). 309 time lapses. (MOV)

Supporting Information
Movie S3 180 degree rotation of 3-D plot shown in Figure 2C. Rotation is done around axis Number of macrophages (MF), while the other two axes are diameter of blood vessels (Dmax) and percent collagen fibers (%Col). Blue balls represent microenvironments that enable fast locomotion and red balls those enabling slow-locomotion. (MOV) Movie S4 Representative example (related to Figure 4A