Extracellular Matrix regulates the morphodynamics of 3D migrating cancer cells

Cell shape is linked to cell function. The significance of cell morphodynamics, namely the temporal fluctuation of cell shape, is much less understood. Here we study the morphodynamics of MDA-MB-231 cells in type I collagen extracellular matrix (ECM). We systematically vary ECM physical properties by tuning collagen concentrations, alignment, and gelation temperatures. We find that morphodynamics of 3D migrating cells are externally controlled by ECM mechanics and internally modulated by Rho-signaling. We employ machine learning to classify cell shape into four different morphological phenotypes, each corresponding to a distinct migration mode. As a result, we map cell morphodynamics at mesoscale into the temporal evolution of morphological phenotypes. We characterize the mesoscale dynamics including occurrence probability, dwell time and transition matrix at varying ECM conditions, which demonstrate the complex phenotype landscape and optimal pathways for phenotype transitions. In light of the mesoscale dynamics, we show that 3D cancer cell motility is a hidden Markov process whereby the step size distributions of cell migration are coupled with simultaneous cell morphodynamics. We also show that morphological phenotype transitions facilitate cancer cells to navigate non-uniform ECM such as traversing the interface between matrices of two distinct microstructures. In conclusion, we demonstrate that 3D migrating cancer cells exhibit rich morphodynamics that is regulated by ECM mechanics, Rho-signaling, and is closely related with cell motility. Our results pave the way to the functional understanding and mechanical programming of cell morphodynamics for normal and malignant cells.


I. INTRODUCTION
Shape defines the cell. In the 1677 book Micrographia, Robert Hooke showed sections within a herbaceous plant under a microscope. The shape of those sections resembles cells in a monastery, so he named the structures cells [1]. Many breakthroughs followed Hooke's discovery, from the cell theory of Schwann and Schleiden, to the theory of tissue formation by Remak, Virchow and Kolliker, and the theory of cellularpathologie by Virchow, all of which are inspired by observations of cell shapes, or morphology in general [2,3].
In our modern view cell shape is determined by cell function [4,5]. A nerve cell has long branched protrusions for communication with other neurons; while the cuboidal shape of epithelial cells allow them to tile the surface of organs. Loss of characteristic shape, on the other hand, is associated with functional abnormality. Thus morphological characterization has been an important tool for diagnosis such as in red blood cell disease [6], neurological disease [7], and cancer [8,9]. More recently, cell shape analysis is boosted by techniques from computer vision. As a result, it becomes possible to obtain high content information of cellular states from morphological data alone [9][10][11][12][13].
While most research focuses on the static cell morphology, the dynamic fluctuation of cell shape is much less understood [14,15]. However, shape fluctuation -namely morphodynamics, is of central importance for dynamic * Correspondence to sunb@onid.orst.edu cellular functions. The abnormal diffusion of small protrusions -microvilli -on the surface of a T cell allows the T cell to efficiently scan antigen-presenting surfaces [16]. For a migrating cancer cell, morphodynamics drives the motility of the cell in many ways similar to our body frame movements that enable swimming. In fact, just as there are different swimming styles, cancer cells have been observed to execute multiple programs -migration modes -during invasion in 3D tissue space [17]. Each mode has distinct signatures of morphology and morphodynamics, and are usually classified based on cell morphology as filopodial, lamellipodial, lobopodial, blebbing, and actin-enriched leading edge [18]. Cancer cell migration modes is controlled by intracellular signaling such as the Rho-ROCK-myosin pathways [19,20], and extracellular factors such as the elasticity, and degradability of the extracellular matrix (ECM) [18,21]. The ability of a cancer cell to switch between migration modes is important for tumor prognosis. Many therapies, such as MMP inhibitors that target a particular mode of cell motility, fail to stop tumor metastasis largely because cells take other available migration programs [22,23].
In this paper, we study the morphodynamics of MDA-MB-231 cells, a highly invasive human breast cancer cell line, in 3D collagen matrices. We devise machine learning techniques to classify cell shapes into morphological phenotypes that correspond to known migration modes. This approach provides a mesoscale mapping of cell morphodynamics into transitions among morphological phenotypes. We find individual cells are capable of rapidly sampling multiple morphological phenotypes, implying spontaneous migration mode transitions. We find ECM mechanics coupled with cell mechanosensing pathways regulate the stability and transition rates between morphological phenotypes. We also find that such transitions facilitate cancer cells to navigate ECM with inherent structural and mechanical heterogeneity. Our results reveal 3D cancer cell migration as a hidden Markov process and morphodynamics contribute to the changes of motility by ECM physical cues. The mean square displacement (σ 2 ) of selected cell geometric measures. Dots: experimental measurements. Solid lines: linear fit. Dashed lines: 95% prediction interval. Here the form factor is defined as perimeter 2 /area. Curl is defined as the ratio between the major axis length and skeletonized contour length. This figure is prepared with Mathlab R2020a (www.mathworks.com) and ImageJ (https://imagej.net).

II. RESULTS
We find 3D migrating cancer cells demonstrate rapid shape fluctuations ( Fig. 1(A-B)). In order to quantify the cell morphodynamics, we take time-lapse fluorescent images of MDA-MB-231 cells migrating in collagen matrices. The GFP-labeled cells typically stay within the focal depth of the objective lens (20X, NA 0.7) for 10-20 hours, while we obtain 2D cell images at a rate of 4 frames per hour (see SI Appendix section S1). After binarization and segmentation, we compute a total of twenty-one geometric measures which collectively quantify the shape of a cell (see SI Appendix section S2). These geometric measures characterize cell size (such as area and perimeter), deviation from circle (such as aspect ratio and form factor), surface topography (such as solidity), and backbone curvature (such as curl -the ratio between the major axis length and skeletonized contour length).
The morphodynamics of a cell manifests itself as a random walk in the geometric shape space concurrent with its motility in the 3D matrix (Fig. 1). However, unlike the real space motility that is slightly superdiffusive [24], we find cell morphodynamics is subdiffusive in the geometric shape space ( Fig. 1C and SI Appendix section S3). The subdiffusivity suggests physical barriers that are present both intrinsic to the cells and from the 3D ECM. Indeed, we find cells moving on 2D surface exhibit faster shape fluctuations than cells embedded in 3D ECM. Still, on flat surfaces cells show subdiffusive random walks in the geometric shape space and superdiffusive walks in real space (SI Appendix section S3).
After quantitatively demonstrating the shape fluctuations of migrating cancer cells, we next investigate cell morphodynamics at a mesoscale that allows us to gain insights on cell migration modes. This is possible because different migration modes are associated with distinct characteristic cell morphologies ( Fig. 2A) [17,18]. Using 3800 manually labeled single cell images, we have trained machine classifiers that classify cell morphology into four morphological phenotypes based on and named after their corresponding migration modes.
We consider four morphological phenotypes including two amoeboidal ones: actin-enriched leading edge (or AE in short) and small blebbing (BB); as well as two mesenchymal ones: filopodial (FP) and lamellipodial (LA). These morphological phenotypes are associated with known migration modes exhibiting characteristic molecular fingerprints. For instance AE and BB cells do not rely on cell-ECM adhesions during migration. For BB cells, their cortical stress continuously drives the formation of rounded blebs at the cell membrane [25][26][27]. AE cells, on the other hand, demonstrate elevated actin polymerization that drive sharp protrusions [28,29]. Neither AE nor BB cells show clear polarization of cytoskeleton and cytoskeleton-associated proteins. In contrast, FP and LA cells exhibit strong cell-ECM adhesions. The filopodial cells consist of distinguishable F-actin bundles extending across the polarized cell body [30,31], while the lamellipodial cells feature fan-shaped leading edges of migration [32,33]. See also Fig. S19 for characteristic Factin subcellular structures. Of note, another migration mode, namely lobopodial or nuclear piston mode, has not been observed in our experiments which is consistent with previous reports [34].
Once the classifier is trained, morphological phenotypes are determined automatically from a cell image if a particular phenotype receives more than 60% probability score (Fig. 2B). For a small fraction of cells (≈ 10%), none of the four phenotypes receive more than 60% probability score, we consider these cells to be in an intermediate state.
We have trained two classifiers (see SI Appendix S4).  The cell images are quantified using a total of 21 geometric measures such as area, solidity, and aspect ratio. With 3800 manually labeled single cell images we have trained a supported vector machine (SVM) to calculate probability scores (Val) for a cell to belong to each morphological phenotypes (classes). We assign a cell to the class C with the maximum score (M ax(C, val)), if this maximum score is greater than a threshold of 60% (M ax(C, val) >0.6). We consider a cell to be at an intermediate state if none of the four classes have a score higher than 0.6. In (B) a sample cell image is classified as a lamellipodial cell (LA), because LA class has a score of greater than 0.6. (C) To better visualize the high dimensional geometric measures, we apply t-SNE method to generate a 3D projection of the geometric cell shape space. The first one is based on support vector machines (SVM [35,36]) involving 21 geometric measures of binarized cell images. The second one is based on Random-Forest model using the same geometric properties [37]. The two classifiers agree with each other well on test data sets (90% overlapping). The SVM classifier particularly has a higher success rate of classifying unseen data (88%). In the following we mainly report the results from SVM algorithm. We also employ the t-SNE algorithm to reduce the dimension of the (21-dimensional) geometric shape space to facilitate visualization of cellular morphodynamics. As shown in Fig. 2C, unseen data (15,000 data points) belonging to different morphological phenotypes form separable clusters in the embedding space. By applying the SVM classifier to time lapse record-ings of 3D migrating MDA-MB-231 cells we find cells spontaneously make transitions among different morphological phenotypes. Fig. 3A shows snapshots of a typical cell. The cell switches directly from blebbing (B) mode to lamellipodial (L) mode via intermediate state (I). Therefore using machine learning technique we map cell morphodynamics into transitions between morphological phenotypes, or their associated migration modes. In order to understand the mechanisms underlying cell morphological phenotype transitions, we examine the effects of manipulating Rho/ROCK-signaling, which is a master regulator that determines the mechanical state of a cell. Rho/ROCK-signaling controls key aspects of cell morphogenesis and migration, such as actomyosin contractility, actin polymerization, cell-cell and cell-ECM adhesion [38]. By regulating the cell mechanotransduction pathways, Rho/ROCK-signaling has also been shown to control 3D cell migration phenotype plasticity [20].
We apply Y27632 [39], a potent Rho kinase (ROCK) inhibitor, and CN03 [40], a Rho-activator to MDA-MB-231 cells cultured in collagen ECM (see SI appendix S5). Y27632 reduces actomyosin contractility, promoting transitions from blebbing to mesenchymal phenotypes [41]. On the other hand, CN03 elevates myosin II activity, leading to retraction of filopodia to rounded cell shapes (Fig. 3B). These results are consistent with previous reports on the molecular control of cell migration modes by Rho/ROCK-signaling [18].
While previous studies focus on the end points of manipulating Rho/ROCK-signaling, morphodynamic analysis offers insights to the transition paths between migration modes. In particular, we take advantage of a modified t-SNE algorithm [42], which projects a cell image in the embedding space defined by the training sets ( Fig. 2C, SI appendix S4). This approach allows us to map the continuous shape change of a cell as a trajectory in the embedded shape space. Similar approaches have also been employed previously in studying complex body movements of other organisms such as fruit flies, where transition paths between different fly behaviors can be visualized [43].
Tracking the mesoscale morphodynamics of MDA-MB-231 cells under pharmacological perturbations, we find up and down regulation of Rho/ROCK-signaling do not lead to a reversal of morphodynamic trajectories. In particular, when treated with Y27632 blebbing cells turn to filopodial or lamellipodial via strongly converging trajectories most of which first visiting AE states (see also SI appendix S5). Fig. 3C shows two representative trajectories. AE state exhibits weak cell-ECM adhesions and F-actin rich protrusions [44,45]. Our results suggest AE states mediate Rho/ROCK-signaling controlled transition from amoeboidal to mesenchymal motility. On the other hand, CN03 treatment causes the majority of mesenchymal cells to switch to blebbing modes. However, without going through AE states, CN03 leads to strongly fluctuating and diverging trajectories observed from multiple cells (Fig. 3D shows two representative trajectories, see also SI appendix S5).
Since Rho/ROCK-signaling is employed by cells to sense ECM physical properties, we next investigate the external control of mesoscale cell morphodynamics. In particular, we focus on the role of ECM physical properties in regulating cell morphological phenotype transitions. In order to control the microstructure of collagen matrices we apply three methods as shown in Fig.  4(A-D) (see also SI appendix S6). First, increasing gelation temperature from room temperature (RT) to 37 • C, while keeping collagen density at [col]= 1.5 mg/mL significantly reduces fiber length and pore size (Fig. 4B). Second, increasing collagen density to [col]= 3.0 mg/mL while keeping gelation temperature at room temperature moderately reduces pore size, preserves a clear fibrous structure, and increases stiffness (Fig. 4C). Finally, keep-ing gelation at room temperature and [col]=1.5 mg/mL, while generating an unidirectional flow field during gelation leads to aligned collagen fibers in the ECM. This method creates strong anisotropy in the ECM.
We find the occurrence probability (population fraction) of different morphological phenotypes are remarkably different at different ECM conditions. As shown in Fig. 4E, increasing gelation temperature does not affect the probability of AE and LA cells. However, the homogeneous matrix microstructure at 37 • C significantly reduces the fraction of FP cells from 43% to 25%, while increases fraction of BB cells from 15% to 25%. Compared with increasing gelation temperature, doubling collagen concentration leads to less dramatic changes of the ECM microstructure. Correspondingly, only moderate changes of phenotype probabilities are observed. On the other hand, when matrix anisotropy is increased by aligning ECM fibers, we find significant shift of cells from amoeboid phenotypes to filopodial mode. Taken together, these results show that ECM heterogeneity and anisotropy determine the probability of different morphological phenotypes.
We have also examined the stability of each morphological phenotype by measuring the average dwell timeduration of a cell to stay continuously in a morphological phenotype before transition to another (SI appendix S7). As shown in Fig. 4F, in all three cases manipulating ECM physical properties moderately increase the dwell time of all four morphological phenotypes. Therefore the changes in the phenotype probability can not be explained by the phenotype stability alone, and in some cases move in opposite trend from the dwell times observed. For instance, occurrence probability of filopodial cells is lowest at 37 • C, which happens to be the condition where filopodial cells show longest dwell time. As we will show later, such discrepancies can be attributed to the detailed transition paths between the phenotypes.
To reveal further details of morphological phenotype dynamics, we have computed the phenotype transition matrix: rates r that characterize the probability of transitions per hour between any two phenotypes ( Fig. 4G-J (see also SI appendix S7). While the rates vary dramatically for different ECM conditions (arrows in Fig.  4G-J), we notice several remarkable common features. First, direct transitions along FP -BB path rarely happen (r <0.03 hr −1 ). Instead, amoeboidal -mesenchymal transitions are primarily mediated by LA and AE states, presumably by turning cell-ECM adhesion on and off. On the other hand, transitions within the amoeboidal (AE -BB) and mesenchymal (FP -LA) modes are frequent, and the rates can go up to 1 per hour. Finally, while the morphological phenotype transitions are intrinsically non-equilibrium processes, probability fluxes between states are generally very small (SI appendix S7). This means that an approximate detailed balance exists among morphological phenotypes. In comparison with other nonequilibium stationary processes at mesoscale [46], we speculate morphological phenotype transitions are not gated by active processes such as ATP consumption. The transition rates also offer insights to understand the ECM-dependence of the fraction of cells in each morphological phenotype (Fig. 4E). For instance, as gelation temperature increases from RT to 37 • C, rates from AE to FP decreases by 52 percent, and rates from BB to AE decreases by 22 percent (Fig. 4G and Fig. 4H). As a result, we observe more blebbing cells and less filopodial cells in collagen matrices prepared at 37 • C. This is consistent with the mechanical mechanism of blebbing formation [47,48]. Blebs form when actomyosin contractility exceeds the binding between cortical actin and cell membrane. A blebbing cell turns to AE when actin polymerization causes sharp protrusion on the membrane. Our results suggest that homogeneous collagen ECM favors BB phenotype to AE, likely due to the reduced protrusive force associated with actin polymerization.
Conversely, as ECM becomes more anisotropic ( Fig.  4G and Fig. 4J), the transition rate from LA to FP increases as much as 27 percent, while rate from AE to BB decrease by 44 percent. Together, these altered rates lead to a significant fraction of blebbing cells turning to filopodial as shown in Fig. 4E. Filopodial protrusions consist of elongated F-actin bundles supported by elevated actin polymerization and cross-linking by Ena/VASP proteins [49]. Our results suggest that the mechanical barrier separating filopodia and blebbing protrusions is too high for actomyosin contractility to overcome directly. Instead, a blebbing cell turning into a filopodial one has to first transform into AE or LA states.
Because the morphological phenotype of a cell is linked to its 3D migration mode, we next investigate if the invasion potential of MDA-MB-231 cells depends on the mesoscale morphodynamics. Due to the short dwell times for each morphological phenotype, we only consider two coarse-grained classes of morphologies: mesenchymal (ME), which consists of FP and LA states; and amoeboidal (AM), which consists of AE and BB states. In particular, we measure for short time scales the step size distributions and for longer time scales the mean square displacement of the cells in randomly aligned collagen matrices gelled at room temperature (Fig. 5).
Interestingly we find the steps are better described by a log-normal, rather than Gaussian distribution (SI appendix S8) due to frequent large steps. Fig. 5A shows the mean and variance of the fitting parameters. It is clear that the steps in physical space are coupled with the corresponding mesoscale dynamics. For cells that dwell in the amoeboidal class, both mean and variance of the steps are the smallest. Correspondingly, the mean square displacement of amoeboid cells have a small slope, corresponding to an effective diffusivity of 6 µm 2 /hour (for each spatial dimension, Fig. 5B). On the other hand, cells make larger steps when dwelling in the mesenchymal class, and the effective diffusivity increases by three-fold to 19 µm 2 /hour. The observed higher motility for cells in the ME state than in AE state is consistent with both migratory measurements in vitro [50] and metastasis measurements in vivo [51]. For instance, using a mouse breast cancer model, it is shown that inhibition of mesenchymal phenotype by NEDD9-depletion significantly (by ≈ 50%) reduces the number of circulating tumor cells.
Our analysis shows that not only it is important to distinguish different morphological phenotypes in studying the motility of cancer cells, but also one may need to take into account of phenotype transitions. We find cell migration steps associated with different class-switching events have distinct statistical distributions (Fig. 5A, also see SI appendix S8 for the motility characteristics of cells involving intermediate states). Without accounting for the class-switching events, the weighted average of mean square displacements from mesenchymal, amoeboidal, and intermediate state cells underestimates the observed cell motility by 20% (Fig. 5B, the weighted average MSD curve corresponds to a diffusivity of approximately 15 µm 2 /hour, as compared with 19 µm 2 /hour for full cell trajectories). We also find cell motility depend on the direction of state transitions. For instance, the variance of steps coupled to ME-to-AM transition is almost twice the value for AM-to-ME transition (Fig. 5A, lower  panel). Since morphological phenotype transitions occur spontaneously at single cell level, our results show that morphodynamics and cell motility are closely coupled to determine the invasive potential of cancer cells.
During metastasis a migrating cancer cell must navigate ECM of distinct mechanical properties. Therefore, we next investigate how cell morphodynamics facilitate cell traverse interfaces and adapt to ECM of distinct mechanics. To this end, we create collagen matrices consist of two integrated layers (SI appendix S9). The RT layer is prepared at room temperature that shows a porous fibrous network, and the 37 • C layer is prepared at 37 • C showing a much more homogeneous structure (Fig.  6A). Without additional cues MDA-MB-231 cells randomly navigate the ECM, occasionally traverse the interface to experience a sudden change of ECM physical properties. Over the course of 24 hours we do not observe durotaxis.
Consistent with the corresponding results in uniform ECM, we find the likelihood of observing a filopodial cell is significantly higher in the ECM layer prepared at room temperature, while for blebbing cells the probability is higher in the gel layer prepared at 37 • C (Fig. 6B). The dwell times of phenotypes, on the other hand, follow the same trend of occurrence probabilities but change rather moderately (Fig. 6C).
The shift of phenotype homeostasis once again can be understood from the phenotype transition matrices. To simplify the discussion in Fig. 6D we plot the top four matrix elements of the transition matrices calculated from cells in each of the two layers. In the RT layer, filopodial cell population is enriched by frequent LA-AE AM-AM AM-ME ME-ME ME-AM 0 exchange and LA to FP transition. As cells cross the interface, LA to FP becomes less likely, and the AE-BB pathway is steered to favor blebbing states.
To further quantify the effect of the interface in modulating cell morphodynamics, we calculate spatial frequencies of dwell events (Fig. 6E) and AE-originating transition events (Fig. 6F). We define a coordinate system where the y-axis is along the interface passing x=0 (Fig.  6A). This allows us to combine data from multiple repeating experiments where cell locations are seeded randomly. After aligning the coordinates, we define the spatial frequency of transition (or dwell) events from state i to state j as R(i → j, x), which satisfies Here P (i → j, x) is the probability density of observing event i → j per unit time (1 hour), and M (x) is the cell density (along x-axis). We use a 1-D Gaussian kernel to estimate P (i → j, x) and M (x) (SI appendix S9). As a result, the spatial frequency R(i → j, x) represents the likelihood of a cell to undergo a specific type of transition over unit time (1 hour) as a function of distance to the interface. The spatial frequency of dwell events clearly show that while FP state is increasingly stable into the RT layer, LA and BB states are more stable in the 37 • C layer. AE state, on the other hand, is most stable at the interface (Fig. 6E). Therefore AE state plays a special role in mediating the cell adaptation across distinct ECM layers. Indeed, we find a gradual shift of favorable AEoriginating transitions as distance to the interface varies. The frequency of AE to LA events, the main amoeboidal to mesenchymal path, peaks in the RT layer. AE to BB events, which is mainly responsible of enriching blebbing cells, has peak frequency in the 37 • C layer. Taken together, we find morphological phenotype transitions and the associated migration mode switching are integral parts of cancer cell invasion and adaptation to complex ECM.

III. DISCUSSION
In this paper, we report the morphodynamics of MDA-MB-231 cells in type I collagen ECM as a model system of metastatic cancer cells migrating in 3D tissue. MDA-MB-231 cells rapidly change their geometry, exhibiting a subdiffusive random walk in the geometric shape space. This occurs simultaneously with their superdiffusive walks in the real space (Fig. 1).
The biological significance of the morphodynamics is further demonstrated by classifying cell shapes into morphological phenotypes corresponding to different migration programs (Fig. 2). This allows us to study cell morphodynamics at the mesoscale, in terms of morphological phenotype transitions. Utilizing machine learning and visualization techniques, we show that cell morphodynamics is regulated by Rho/ROCK-signaling (Fig. 3), which is a molecular control hub of cell mechanosensing and force generation. It has been shown previously that Rho/Rac signaling regulates the shift between mesenchymal and amoeboidal motility [19,21]. Our analysis further reveals that instead of favoring a particular mode of motility, perturbations of Rho/ROCK-signaling alter the migration mode transition rates. In particular, down regulating Rho leads to overall amoeboidal-tomesenchymal transition that routes through AE and LA states. Activation of Rho, on the other hand, leads to strongly fluctuating morphodynamics that enriches blebbing cells. The irreversibility of up and down regulating Rho/ROCK-signaling results suggest a complex phenotype landscape that controls 3D cancer cell motility.
We study morphological phenotype transitions in ECM of distinct physical properties and find ECM microstructure modulates the probabilities, dwell times, and transition rates of morphological phenotypes. Collagen matrices with homogeneous structure, as those prepared at higher temperature, enrich the population of blebbing cells. By comparing the transition matrices, we find the enrichment of blebbing cells is directly related with the reduced transition rate from BB to AE state, and also indirectly contributed by the mesenchymal-to-amoeboidal transition through LA and AE states. Similarly, collagen matrices with structural anisotropy enrich the population of filopodial cells. The enrichment is directly attributed to an increased LA to FP rate, and indirectly contributed by the amoeboidal-to-mesenchymal transition mediated by LA and AE states. These results show that it is possible to execute external control of cell morphodynamics (and the corresponding 3D migration modes) through ECM mechanics. Importantly, taking into account of the phenotype transitions allows us to better predict the outcome of manipulating cell migration mode through ECM physical properties [17,52].
In light of the rapid phenotype transitions exhibited by individual cells, 3D cancer cell motility may be considered as a hidden Markov process where each phenotype is associated with characteristic step size distributions (Fig.  5). Specifically, we find steps that occur simultaneously with a phenotype transition have distinct sizes compared with steps that occur while cells dwell in a particular morphological phenotype. This makes morphodynamics a crucial factor in determining the invasive potential of cancer cells. To our knowledge, this aspect has been so far largely overlooked in the literature.
In the lens of a hidden Markov process, morphodynamics may facilitate cancer invasion because phenotype transitions allow cancer cells to search for and commit to a more effective migration program [53]. Using a ECM model consisting of two mechanically distinct layers, we show the cells gradually adjust their morphodynamics as they approach and cross the layer interface (Fig. 6). Therefore morphological phenotype transitions may be essential in cancer cell metastasis by enabling the cells to navigate non-uniform ECM.
The connection between morphological phenotypes and cell migration modes may be further strengthened by incorporating key molecular fingerprints. Using a small data set, we find the machine-classified morphological phenotypes have F-actin structures that are expected from the corresponding migration modes (Fig. S19, [31]). We anticipate that training with a multichannel data set which includes not only cell shape but also molecular profiles such as actin, myosin, and integrin, will improve accuracy in identifying cell migration modes. However, doing so could reduce the throughput and limit the applicability in general imaging settings.
In summary, we demonstrate the morphodynamics of 3D migrating cancer cells as a powerful tool to inspect the internal state and microenvironment of the cells. Investigated at mesoscale, the morphodynamics imply that 3D cancer cell migration is inherently plastic [18]. The plasticity is controlled by the mode transition matrices, rather than a deterministic decision tree [17]. In order to further exploit the information provided by the cell shape fluctuations, future research is needed to decode morphodynamics as a rich body language of cells, and to control morphodynamics as a route of mechanical programming of cell phenotype.

IV. MATERIALS AND METHODS
GFP-labeled MDA-MB-231 cells are purchased from GenTarget Inc and are handled following the vendor's recommendations. High concentration rat tail Type I collagen solutions are purchased from Corning Inc, and the collagen gels are prepared following standard protocol. Cells and collagen fibers are imaged using Leica SPE confocal microscope with fluorescent and reflection channels simultaneously. A streamline of utilizing Im-ageJ , Matlab, and Python scripts are developed to analyze raw images. See SI Appendix S1-S10 for details of 3D cell culture, microscopy, pharmacological treatments, and data analysis.