Extracellular Matrix regulates the morphodynamics of 3D migrating cancer cells
EExtracellular Matrix regulates the morphodynamics of 3D migrating cancer cells
Christopher Z. Eddy, Helena Raposo, Ryan Wong, Fuxin Li, and Bo Sun Department of Physics, Oregon State University, Corvallis OR, 97331 Department of Chemical, Biological, and Environmental Engineering, Oregon State University, Corvallis OR, 97331 School of Electrical Engineering and Computer Science, Oregon State University, Corvallis OR, 97331
Cell shape is linked to cell function. The significance of cell morphodynamics, namely the tem-poral fluctuation of cell shape, is much less understood. Here we study the morphodynamics ofMDA-MB-231 cells in type I collagen extracellular matrix (ECM). We systematically vary ECMphysical properties by tuning collagen concentrations, alignment, and gelation temperatures. Wefind that morphodynamics of 3D migrating cells are externally controlled by ECM mechanics andinternally modulated by Rho-signaling. We employ machine learning to classify cell shape into fourdifferent 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 transitionmatrix at varying ECM conditions, which demonstrate the complex phenotype landscape and op-timal pathways for phenotype transitions. In light of the mesoscale dynamics, we show that 3Dcancer cell motility is a hidden Markov process whereby the step size distributions of cell migrationare coupled with simultaneous cell morphodynamics. We also show that morphological phenotypetransitions facilitate cancer cells to navigate non-uniform ECM such as traversing the interface be-tween matrices of two distinct microstructures. In conclusion, we demonstrate that 3D migratingcancer cells exhibit rich morphodynamics that is regulated by ECM mechanics, Rho-signaling, andis closely related with cell motility. Our results pave the way to the functional understanding andmechanical programming of cell morphodynamics for normal and malignant cells.
INTRODUCTION
Shape defines the cell. In the 1677 book
Micrographia ,Robert Hooke showed sections within a herbaceous plantunder a microscope. The shape of those sections resem-bles cells in a monastery, so he named the structurescells [1]. Many breakthroughs followed Hooke’s discov-ery, from the cell theory of Schwann and Schleiden, to thetheory of tissue formation by Remak, Virchow and Kol-liker, and the theory of cellularpathologie by Virchow, allof which are inspired by observations of cell shapes, ormorphology in general [2, 3].In our modern view cell shape is determined by cellfunction [4, 5]. A nerve cell has long branched protru-sions for communication with other neurons; while thecuboidal shape of epithelial cells allow them to tile thesurface of organs. Loss of characteristic shape, on theother hand, is associated with functional abnormality.Thus morphological characterization has been an impor-tant 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 com-puter vision. As a result, it becomes possible to obtainhigh content information of cellular states from morpho-logical data alone [9–12].While most research focuses on the static cell morphol-ogy, the dynamic fluctuation of cell shape is much less un-derstood. However, shape fluctuation – namely morpho-dynamics, is of central importance for dynamic cellularfunctions. The abnormal diffusion of small protrusions -microvilli - on the surface of a T cell allows the T cellto efficiently scan antigen-presenting surfaces [13]. For a migrating cancer cell, morphodynamics drives the motil-ity of the cell in many ways similar to our body framemovements that enable swimming. In fact, just as thereare different swimming styles, cancer cells have been ob-served to execute multiple programs – migration modes– during invasion in 3D tissue space [14]. Each mode hasdistinct signatures of morphology and morphodynamics,and are usually classified based on cell morphology asfilopodial, lamellipodial, lobopodial, blebbing, and actin-enriched leading edge [15]. Cancer cell migration modesis controlled by intracellular signaling such as the Rho-Rock-myosin pathways [16, 17], and extracellular factorssuch as the elasticity, and degradability of the extracellu-lar matrix (ECM) [15, 18]. The ability of a cancer cell toswitch between migration modes is important for tumorprognosis. Many therapies, such as MMP inhibitors thattarget a particular mode of cell motility, fail to stop tu-mor metastasis largely because cells take other availablemigration programs [19, 20].In this paper, we study the morphodynamics of MDA-MB-231 cells, a highly invasive human breast cancer cellline in 3D collagen matrices. We devise machine learn-ing techniques to classify cell shapes into morphologicalphenotypes that correspond to known migration modes.This approach provides a mesoscale mapping of cell mor-phodynamics into transitions among morphological phe-notypes. We find individual cells are capable of rapidlysampling multiple morphological phenotypes, implyingspontaneous migration mode transitions. We find ECMmechanics coupled with cell mechanosensing pathwaysregulate the stability and transition rates between mor-phological phenotypes. We also find that such transi- a r X i v : . [ q - b i o . CB ] O c t tions facilitate cancer cells to navigate ECM with inher-ent structural and mechanical heterogeneity. p e r i m e t e r ( 𝜇 𝑚 ) a s p ec t r a ti o time (hour)0 25 A B time delay (hour) 𝜎 𝐴 𝑟 𝑒 𝑎 ( 𝜇 𝑚 ) 𝜎 𝑓 𝑜 𝑟 𝑚 𝑓 𝑎 𝑐 𝑡 𝑜 𝑟 ( × − ) 𝜎 𝑐 𝑢 𝑟 𝑙 ( × − ) 𝑦~𝑥 𝑛 𝑛 = 0.6 𝑦~𝑥 𝑛 𝑛 = 0.4 𝑦~𝑥 𝑛 𝑛 = 0.25 C FIG. 1. Three-dimensional migration of MDA-MB-231 cellsis accompanied with significant cell shape fluctuation. (A) Atypical time lapse recording of 25 hours is projected onto asingle image with colors representing time. (B) The real space(x-y plane), and shape fluctuations of 3 cells shown in (A). (C)The mean square displacement ( σ ) of selected cell geometricmeasures. Dots: experimental measurements. Solid lines:linear fit. Dashed lines: 95% prediction interval. Here theform factor is defined as perimeter / area. Curl is definedas the ratio between the major axis length and skeletonizedcontour length. RESULTS
We find 3D migrating cancer cells demonstrate rapidshape fluctuations (Fig. 1(A-B)). In order to quantifythe cell morphodynamics, we take time-lapse fluorescentimages of MDA-MB-231 cells migrating in collagen ma-trices. The GFP-labeled cells typically stay within thefocal depth of the objective lens (20X, NA 0.7) for 10-20hours, while we obtain 2D cell images at a rate of 4 framesper hour (see
SI Appendix section S1). After binarizationand segmentation, we compute a total of twenty-one geo-metric measures which collectively quantify the shape ofa cell (see
SI Appendix section S2). These geometric mea-sures characterize cell size (such as area and perimeter),deviation from circle (such as aspect ratio and form fac-tor), surface topography (such as solidity), and backbonecurvature (such as curl – the ratio between the major axislength and skeletonized contour length).The morphodynamics of a cell manifests itself as a ran-dom walk in the geometric shape space concurrent with its motility in the 3D matrix (Fig. 1). However, unlikethe real space motility that is slightly diffusive, we findcell morphodynamics is subdiffusive in the shape space(Fig. 1C and
SI Appendix section S3). The subdiffusivitysuggests physical barriers that are present both intrinsicto the cells and from the 3D ECM. Indeed, we find cellsmoving on 2D surface exhibit faster shape fluctuationsthan cells embedded in 3D ECM. Still, on flat surfacescells show subdiffusive random walks in shape space andsuperdiffusive walks in real space (SI Appendix sectionS3). V V2 Train Data V1 geometric measures t-SNE embedding Area = 693 µm Solidity = 0.6
Aspect Ratio = 1.1(18 other measures)
SVM Classifier C AE BBFPLA
Val actin-enriched leading edge (AE) blebbing (BB) filopodia (FP) lamellipodia (LA) visualization of training sets
AB C
AEBBFPLA
FIG. 2. Develop a supervised machine learning model to clas-sify cells into morphological phenotypes corresponding to dif-ferent migration modes. (A) MDA-MB-231 cells in 3D colla-gen matrices exhibit multiple morphological phenotypes thatare characteristic of four distinct migration modes: actin-enriched leading edge (AE), small blebbing (BB), filopodial(FP), and lamellipodial (LA). Scale bars are 20 µ m. (B) Thecell images are quantified using a total of 21 geometric mea-sures such as area, solidity, and aspect ratio. With 3800manually labeled single cell images we have trained a sup-ported 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 maximumscore ( Max ( C, val )), if this maximum score is greater than athreshold of 60% (
Max ( C, val ) > After quantitatively demonstrating the shape fluctua-tions of migrating cancer cells, we next investigate cellmorphodynamics at a mesoscale that allows us to gaininsights on cell migration modes. This is possible be-cause different migration modes are associated with dis-tinct characteristic cell morphologies (Fig. 2A) [14, 15].Using 3800 manually labeled single cell images, we havetrained machine classifiers that classify cell morphologyinto four morphological phenotypes based on and namedafter their corresponding migration modes.We consider four morphological phenotypes includingtwo amoeboidal ones: actin-enriched leading edge (or AEin short) and small blebbing (BB); as well as two mes-enchymal ones: filopodial (FP) and lamellipodial (LA).Of note, another migration mode, namely lobopodial ornuclear piston mode, has not been observed in our ex-periments which is consistent with previous reports [21].Once the classifier is trained, morphological phenotypesare determined automatically from a cell image if a par-ticular phenotype receives more than 60% probabilityscore (Fig. 2B). For a small fraction of cells ( ≈ SI Appendix
S4).The first one is based on support vector machines (SVM[22, 23]) involving 21 geometric measures of binarizedcell images. The second one is based on Random-Forestmodel using the same geometric properties [24]. The twoclassifiers agree with each other well on test data sets(90% overlapping). The SVM classifier particularly hasa higher success rate of classifying unseen data (92%).In the following we mainly report the results from SVMalgorithm. The validity of the SVM classifier is also ev-ident from the color-coded t-SNE embedding of unseendata (Fig. 2C, 25,000 data points), as data points belongto different classes (morphological phenotypes) form sep-arable clusters in the embedding space.By applying the SVM classifier to time lapse recordingsof 3D migrating MDA-MB-231 cells we find cells sponta-neously make transitions among different morphologicalphenotypes. Fig. 3A shows snapshots of a typical cell.The cell switches directly from filopodial (F) to lamel-lipodial (L) mode, then to small blebbing (B) mode viaintermediate state (I). Therefore using machine learningtechnique we map cell morphodynamics into transitionsbetween morphological phenotypes, or their associatedmigration modes.In order to understand the mechanisms underlying cellmorphological phenotype transitions, we examine the ef-fects of manipulating Rho-signaling, which is a masterregulator that determines the mechanical state of a cell.To this end, we apply Y27632 [25], a Rho-inhibitor; andCN03 [26], a Rho-activator to MDA-MB-231 cells cul-tured in collagen ECM (see
SI appendix
S5). Y27632reduces actomyosin contractility, promoting transitionsfrom blebbing to mesenchymal phenotypes [27]. On theother hand, CN03 elevates myosin II activity, leadingto retraction of filopodia to rounded cell shapes (Fig.3B). These results are consistent with previous reportson the molecular control of cell migration modes by Rho-
Rho ↑CN03Rho ↓
Y27632Rho ↓Y27632
Rho ↑
CN03 V1 V2 V3 V1 V2V3 A C D
AEBBFP LA B time FIG. 3. Rho-signaling internally controls mesoscale morpho-dynamics of 3D cultured MDA-MB-231 cells. (A) A sampletime series of morphological phenotype. Insets: three snap-shots showing the GFP-labeled cell morphology. Abbrevia-tions: F – filopodial, L – lamellipodial, I – intermediate state,B: small blebbing. (B) Representative morphological changesunder treatment of CN03 and Y27632, which upregulates anddownregulates Rho A respectively. (C-D) Characteristic mor-phodynamic trajectories of cells in the t-SNE embedded shapespace. The trajectories start immediately after introducingY27632 or CN03, and ends after 4 hours of incubating withthe drugs. The forward time directions are shown as thickcurves with arrows as guide to the eyes. Two representa-tive trajectories (one with circular symbols, and another onewith triangular symbols) per each treatment are shown as col-ored symbols connected by black lines, where color representsthe instantaneous phenotype. Unconnected light-colored dotsshow training sets which are the same as in Fig. 2C. Note tobetter visualize the 3D trajectories coordinates have been ro-tated with respect to Fig. 2C. signaling [15].While previous studies focus on the end points of ma-nipulating Rho-signaling, morphodynamic analysis offersinsights to the transition paths between migration modes.In particular, we take advantage of a modified t-SNEalgorithm [28], which projects a cell image in the em-bedding space defined by the training sets (Fig. 2C,
SI appendix
S4). This approach allows us to map thecontinuous shape change of a cell as a trajectory in the(dimensional-reduced) mesoscale morphodynamic space.Similar approaches have also been employed previouslyin studying complex body movements of other organismssuch as fruit flies, where transition paths between differ-ent fly behaviors can be visualized [29].Tracking the mesoscale morphodynamics of MDA-MB-231 cells under pharmacological perturbations, we findup and down regulation of Rho signaling do not leadto reversal morphodynamic trajectories. In particular,when treated with Y27632 blebbing cells turn to filopo-dial or lamellipodial via strongly converging trajectoriesmost of which first visiting AE states (see also
SI ap-pendix
S5). Fig. 3C shows two representative trajecto-ries. AE state exhibits weak cell-ECM adhesions andF-actin rich protrusions [30, 31]. Our results suggestAE states mediate Rho-signaling controlled transitionfrom amoeboidal to mesenchymal motility. On the otherhand, CN03 treatment causes the majority of mesenchy-mal cells to switch to blebbing modes. However, withoutgoing through AE states, CN03 leads to strongly fluctu-ating and diverging trajectories from multiple cells (Fig.3D shows two representative trajectories, see also
SI ap-pendix
S5).We next investigate the external control of cell mor-phodynamics. In particular, we focus on the role ofECM physical properties in regulating cell morphologi-cal phenotype transitions. In order to control the mi-crostructure of collagen matrices we employ three meth-ods as shown in Fig. 4(A-D) (see also
SI appendix
S6).First, increasing gelation temperature from room tem-perature (RT) to 37 ◦ C, while keeping collagen densityat [col]= 1.5 mg/mL significantly reduce fiber length andpore size (Fig. 4B). Second, increasing collagen densityto [col]= 3.0 mg/mL while keeping gelation temperatureat room temperature moderately reduces pore size, pre-serves clear fibrous structure, and increases stiffness (Fig.4C). Finally, keeping gelation at room temperature and[col]=1.5 mg/mL, while generating an unidirectional flowfield during gelation leads to aligned collagen fibers inthe ECM. This method creates strong anisotropy in theECM.We find the occurrence probability (or population frac-tion) of different morphological phenotypes are remark-ably different at different ECM conditions. As shown inFig. 4E, increasing gelation temperature does not affectthe probability of AE and LA cells. However, the ho-mogeneous matrix microstructure at 37 ◦ C significantlyreduces the fraction of FP cells from 43% to 25%, whileincreases fraction of BB cells from 15% to 25%. Com-pared with increasing gelation temperature, doubling col-lagen concentration leads to less dramatic changes ofthe ECM microstructure. Correspondingly, only mod-erate changes of phenotype probabilities are observed.On the other hand, when matrix anisotropy is increasedby aligning ECM fibers, we find significant shift of cellsfrom amoeboid phenotypes to filopodial mode. Taken to-gether, these results show that ECM heterogeneity andanisotropy determine the probability of different morpho-logical phenotypes.We have also examined the stability of each morpho-logical phenotypes by measuring the average dwell time– duration of a cell to stay continuously in a morpholog-ical phenotype before transition to another (
SI appendix
S7). As shown in Fig. 4F, in all three cases manipu-lating ECM physical properties moderately increase the dwell time of all four morphological phenotypes. There-fore the changes in the phenotype probability can not beexplained by the phenotype stability alone, and in somecases move in opposite trend from the dwell times ob-served.To reveal further details of morphological phenotypedynamics, we have computed the phenotype transitionmatrix: rates r that characterize the probability of tran-sitions per hour between any two phenotypes (Fig. 4G-J(see also SI appendix
S7). While the rates vary dra-matically for different ECM conditions (arrows in Fig.4G-J), we notice several remarkable common features.First, direct transitions along FP - BB path rarely hap-pens ( r < − ). Instead, amoeboidal - mesenchymaltransitions 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 fre-quent, and the rates can go up to 1 per hour. Finally,while the morphological phenotype transitions are intrin-sically non-equilibrium processes, probability fluxes be-tween states are generally very small ( SI appendix
S7).This means that an approximate detailed balance existsamong morphological phenotypes. In comparison withother nonequilibium stationary processes at mesoscale[32], we speculate morphological phenotype transitionsare not gated by active processes such as ATP consump-tion.The transition rates also offer insights to understandthe ECM-dependence of the fraction of cells in each mor-phological phenotype (Fig. 4E). For instance, as gelationtemperature increases from RT to 37 ◦ C, rates from AEto FP decreases by 52 percent, and rates from BB toAE decreases by 22 percent (Fig. 4G and Fig. 4I). Asa result, we observe more blebbing cells and less filopo-dial cells in collagen matrices prepared at 37 ◦ C. Thisis consistent with the mechanical mechanism of blebbingformation [33, 34]. Blebs form when actomyosin contrac-tility exceeds the binding between cortical actin and cellmembrane. A blebbing cell turns to AE when actin poly-merization causes sharp protrusion on the membrane.Our results suggest that as collagen ECM loses struc-ture heterogeneity, actin polymerziation is less effectiveto drive the transition from BB to AE states.Conversely, as ECM becomes more anisotropic (Fig.4G and Fig. 4J), the transition rate from LA to FP in-creases as much as 27 percent, while rate from AE to BBdecrease by 44 percent. Together, these altered rates leadto a significant fraction of blebbing cells turning to filopo-dial as shown in Fig. 4E. Filopodial protrusions consistof elongated F-actin bundles supported by elevated actinpolymerization and cross-linking by Ena/VASP proteins[35]. Our results suggest that the mechanical barrier sep-arating filopodia and blebbing protrusions is too high foractomyosin contractility to overcome directly. Instead,a blebbing cell turning into a filopodial one has to first
AE FP BB LA0123 Dw e ll T i m e ( hou r) AE FP BB LA0123 Dw e ll T i m e ( hou r) AE FP BB LA0123 Dw e ll T i m e ( hou r) AE FP BB LA0123 Dw e ll T i m e ( hou r) AE FP BB LA0123 Dw e ll T i m e ( hou r) AE FP BB LA00.20.40.6 O b s e r v a t i on P r obab ili t y AB E G IH J rate ( ℎ𝑟 −1 ) 10 random, RT AE FP BB LA00.20.40.6 O b s e r v a t i on P r obab ili t y AE FP BB LA00.20.40.6 O b s e r v a t i on P r obab ili t y AE FP BB LA00.20.40.6 O b s e r v a t i on P r obab ili t y AE FP BB LA00.20.40.6 O b s e r v a t i on P r obab ili t y AE FP BB LA00.20.40.6 O b s e r v a t i on P r obab ili t y AE FP BB LA00.20.40.6 O b s e r v a t i on P r obab ili t y AE FP BB LA00.20.40.6 O b s e r v a t i on P r obab ili t y AE FP BB LA00.20.40.6 O b s e r v a t i on P r obab ili t y AE FP BB LA00.20.40.6 O b s e r v a t i on P r obab ili t y AE FP BB LA00.20.40.6 O b s e r v a t i on P r obab ili t y AE FP BB LA00.20.40.6 O b s e r v a t i on P r obab ili t y AE FP BB LA00.20.40.6 O b s e r v a t i on P r obab ili t y AE FP BB LA p r ob a b ilit y F d w e ll ti m e ( h r) collagen densityx2RT to 37°C CD AE FP BB LA AE FP BB LA
FIG. 4. Physical properties of collagen ECM regulate the morphological phenotype homeostasis of 3D migrating MDA-MB-231 cells. (A-D) Confocal reflection images and pseudo colored MDA-MB-231 cells for collagen matrices prepared at varyingconditions. Scale bars: 20 µ m. A: collagen ECM prepared at room temperature (RT, or 25 ◦ C) and collagen concentration of[ col ] = 1.5 mg/mL. B: collagen ECM prepared at 37 ◦ C and [col] = 1.5 mg/mL. C: collagen ECM prepared at RT and [ col ] =3.0 mg/mL. D: collagen ECM prepared with flow-aligned collagen fibers. (E) Fraction of cells in each morphological phenotype.8,000 single cell images are analyzed under each ECM condition. (F) Dwell time of cells in each morphological phenotype.Errorbars in (E-F) represent 95% confidence intervals calculated from 1000 bootstrap iterations. (G-J): The transition matrix– morphological phenotype transition rates under varying ECM conditions. G: collagen ECM prepared at room temperatureand [ col ] = 1.5 mg/mL. H: collagen ECM prepared at 37 ◦ C and [ col ] = 1.5 mg/mL. I: collagen ECM prepared at RT and [col]= 3.0 mg/mL. J: collagen ECM prepared with flow-aligned collagen fibers. Under each ECM condition a total of more than2,000 hours of single cell trajectories are analyzed. transform into AE or LA states.Because the morphological phenotype of a cell is linkedto its 3D migration mode, we next investigate if the in-vasion potential of MDA-MB-231 cells depends on themesoscale morphodynamics. Due to the short dwelltimes for each morphological phenotype, we only considertwo coarse-grained classes of morphologies: mesenchymal(ME), which consists of FP and LA states; and amoe-boidal (AM), which consists of AE and BB states. Inparticular, we measure for short time scales the step sizedistributions and for longer time scales the mean squaredisplacement of the cells in randomly aligned collagenmatrices geled at room temperature (Fig. 5).Interestingly we find the steps are better described bya log-normal, rather than Gaussian distribution (
SI ap-pendix
S8) due to frequent large steps. Fig. 5A showthe mean and variance of the fitting parameters. It isclear that the steps in physical space are coupled withthe corresponding mesoscale dynamics. For cells thatdwell in the amoeboidal class, both mean and varianceof the steps are the smallest. Correspondingly, the meansquare displacement of amoeboid cells have a small slope,corresponding to an effective diffusivity of 8 µ m /hour(for each spatial dimension, Fig. 5B). On the otherhand, cells make larger steps when dwell in the mes-enchymal class, and the effective diffusivity increases bythree-fold to 26 µ m /hour. Importantly, when cell mi-gration is simultaneously coupled with phenotype tran-sitions, the class-switching steps have distinct statisticaldistributions (Fig. 5A). Our analysis shows that not only it is important todistinguish different morphological phenotypes in study-ing the motility of cancer cells, but also one may needto take into account of phenotype transitions. For in-stance, without accounting for the class-switching events,the weighted average of mean square displacements frommesenchymal and amoeboidal cells underestimates theobserved cell motility by 15% (Fig. 5B, the weighted av-erage MSD curve corresponds to a diffusivity of approxi-mately 20 µ m /hour, as compared with 24 µ m /hour forfull cell trajectories). Since phenotype transitions occurrapidly at single cell level regardless of ECM concentra-tion, porosity and rigidity, we conclude that mesoscaledynamics contributes to determine the invasive potentialof cancer cells.During metastasis a migrating cancer cell must navi-gate ECM of distinct mechanical properties. Thereforewe next investigate how cell morphodynamics facilitatecell traverse interfaces and adapt to ECM of distinct me-chanics. To this end, we create collagen matrices consistof two integrated layers ( SI appendix
S9). The RT layeris prepared at room temperature that shows a porousfibrous 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 ran-domly navigate the ECM, occasionally traverse the in-terface to experience a sudden change of ECM physicalproperties. Over the course of 24 hours we do not observedurotaxis.Consistent with the corresponding results in uniform A M - A M A M - M E M E - A M M E - M E m ( m ) A M - A M A M - M E M E - A M M E - M E ( m ) 𝑚 ( µ 𝑚 ) 𝜎 ( µ 𝑚 ) lag time (hours)0.5 1 1.5 2 𝜎 𝑥 𝑦 ( µ 𝑚 ) A B
AMMEaveragefull trajectory
FIG. 5. Cancer cell migration in 3D ECM is coupled withmorphological phenotype and phenotype transitions. (A & B)Means ( µ ) and variances ( σ ) of the step size distributions byfitting the experimental measurements with log-normal dis-tribution. The steps are categorized based on morphologicalphenotype (coarse-grained as mesenchymal and amoeboidal)dynamics. ME: abbreviation for mesenchymal, which includesFP and LA states. AM: abbreviation for amoeboidal, whichincludes BB and AE states. If a cell makes a step by start-ing from AM state to ME state, the step is categorized as anAM-ME step. The starting and ending frames of steps areseparated by 15 minutes. In (A-B) errorbars show the 95%confidence interval of fitted parameters. (C) Real-space (2Dprojection) mean square displacements (MSD) of cells. AM:the MSD of cells dwelling in the AM state. ME: the MSDof cells dwelling in the ME state. average: the weighted av-erage of mean square displacements of AM and ME dwellingcells. The average is based on the occurrence fraction of AMand ME cells. Full trajectory: the MSD obtained from entirecell trajectories regardless of the morphological phenotypes.Shaded areas in C show SEM. In (A-C) a total of 1974 hoursof single cell trajectories are analyzed. The ECM are pre-pared at room-temperature with a concentration [col] = 1.5mg/mL. ECM, we find the likelihood of observing a filopodial cellis significantly higher in the ECM layer prepared at roomtemperature, while for blebbing cells the probability ishigher in the gel layer prepared at 37 ◦ C (Fig. 6B). Thedwell times of phenotypes, on the other hand, follow thesame trend of occurrence probabilities but change rathermoderately (Fig. 6C).The shift of phenotype homeostasis once again canbe understood from the phenotype transition matrices.To simplify the discussion in Fig. 6D we plot the topfour matrix elements of the transition matrices calcu-lated from cells in each of the two layers. In the RT layer,filopodial cell population is enriched by frequent LA-AEexchange and LA to FP transition. As cells cross theinterface, LA to FP becomes less likely, and the AE-BBpathway is steered to favor blebbing states.To further quantify the effect of the interface in modu-lating cell morphodynamics, we calculate spatial frequen-cies of dwell events (Fig. 6E) and AE-originating tran-sition events (Fig. 6F). We define a coordinate systemwhere the y-axis is along the interface passing x =0 (Fig. 6A). This allows us to combine data from multiple repeat-ing experiments where cell locations are seeded randomly.After aligning the coordinates, we define the spatial fre-quency of transition (or dwell) events from state i to state j as R ( i → j, x ), which satisfies P ( i → j, x ) = R ( i → j, x ) M ( x ) (1)Here P ( i → j, x ) is the probability density of observingevent i → j per unit time (1 hour), and M ( x ) is the celldensity (along x -axis). We use a 1-D Gaussian kernel toestimate P ( i → j, x ) and M ( x ) ( SI appendix
S9). As aresult, the spatial frequency R ( i → j, x ) represents thelikelihood of a cell to undergo a specific type of transitionover unit time (1 hour) as a function of distance to theinterface.The spatial frequency of dwell events clearly show thatwhile 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 in-terface (Fig. 6E). Therefore AE state plays a specialrole in mediating the cell adaptation across distinct ECMlayers. Indeed, we find a gradual shift of favorable AE-originating transitions as distance to the interface varies.The frequency of AE to LA events, the main amoeboidalto mesenchymal path, peaks in the RT layer. AE toBB events, which is mainly responsible of enriching bleb-bing cells, has peak frequency in the 37 ◦ C layer. Takentogether, we find morphological phenotype transitionsand the associated migration mode switching are integralparts of cancer cell invasion and adaptation to complexECM.
DISCUSSION
In this paper, we report the morphodynamics of MDA-MB-231 cells in type I collagen ECM as a model systemof metastatic cancer cells migrating in 3D tissue. MDA-MB-231 cells rapidly change their geometry, exhibiting asubdiffusive random walk in the shape space. This occurssimultaneously with their superdiffusive walks in the realspace (Fig. 1).The biological significance of the morphodynamics isfurther demonstrated by classifying cell shapes into mor-phological phenotypes corresponding to different migra-tion programs (Fig. 2). This allows us to study cellmorphodynamics at the mesoscale, in terms of morpho-logical phenotype transitions. Utilizing machine learningand visualization techniques, we show that cell morpho-dynamics is regulated by Rho-signaling (Fig. 3), which isa molecular control hub of cell mechanosensing and forcegeneration. It has been shown previously that Rho/Racsignaling regulates the shift between mesenchymal andamoeboidal motility [16, 18]. Our analysis further re-veals that instead of favoring a particular mode of motil-ity, perturbations of Rho signaling alter the migration p r ob a b ilit y 𝜇𝜇𝑚𝑚 AE→AEBB→BBLA→LAFP→FP 𝑇𝑇 𝑑𝑑 𝑑𝑑 𝑑𝑑 𝑑𝑑𝑑𝑑 ( h r) r a t e ( ℎ 𝑟𝑟 − ) BC D EF
LABBFPAE
AEAEBBBBFPFPLALA 𝑅𝑅 𝑖𝑖 → 𝑖𝑖 ( ℎ 𝑟𝑟 − ) 𝑅𝑅 𝐴𝐴 𝐴𝐴 → 𝑋𝑋 ( ℎ 𝑟𝑟 − ) 𝜇𝜇𝑚𝑚 A BBLAAE
AERT37°C 𝑥𝑥𝑦𝑦 𝑥𝑥𝑥𝑥
FIG. 6. Morphological phenotype transition facilitates cell migration in heterogeneous ECM. (A) Time-lapse projection of3D migrating MDA-MB-231 cells navigating engineered heterogeneous ECM. The ECM contains two adjacent layers that areprepared at room temperature (RT) and 37 ◦ C respectively. A confocal reflection image shows the ECM structure next to theinterface (dashed line, y -axis). Scale bar: 50 µ m. (B) Fraction of cells of each morphological phenotype in both sides of theinterface. (C) Dwell time of cells of each morphological phenotype in both sides of the interface. (D) Phenotype transitionrates in both sides of the interface. (E) Spatial frequency of of dwell events. (F) Spatial frequency of AE to AE, AE to LA andAE to BB events. See main text for the definition of spatial frequency R ( i → j, x ). In (E-F) dashed lines indicate the interface( x =0) separating the 37 ◦ C gel (left) and RT gel (right). A total of 3,800 hours of single cell recordings from three independentexperiments have been used to calculate the results in (B-F). mode transition rates. In particular, down regulatingRho leads to overall amoeboidal-to-mesenchymal transi-tion that routes through AE and LA states. Up regula-tion of Rho, on the other hand, leads to strongly fluctu-ating morphodynamics that enriches blebbing cells. Theirreversibility of up and down regulating Rho signalingresults suggest a complex phenotype landscape that con-trols 3D cancer cell motility.We study morphological phenotype transitions in ECMof distinct physical properties and find ECM microstruc-ture modulates the probabilities, dwell times, and tran-sition rates of morphological phenotypes. Collagen ma-trices with homogeneous structure, as those prepared athigher temperature, enrich the population of blebbingcells. By comparing the transition matrices, we find theenrichment of blebbing cells is directly related with thereduced transition rate from BB to AE state, and also in-directly contributed by the mesenchymal-to-amoeboidaltransition through LA and AE states. Similarly, collagenmatrices with structural anisotropy enrich the populationof filopodial cells. The enrichment is directly attributedto an increased LA to FP rate, and indirectly contributedby the amoeboidal-to-mesenchymal transition mediatedby LA and AE states. These results show that it is pos-sible to execute external control of cell morphodynam-ics (and the corresponding 3D migration modes) throughECM mechanics. Importantly, taking into account of thephenotype transitions allows us to better predict the out-come of manipulating cell migration mode through ECMphysical properties [14, 36].In light of the rapid phenotype transitions exhibited byindividual cells, 3D cancer cell motility may be consid- ered as a hidden Markov process where each phenotype isassociated with characteristic step size distributions (Fig.5). Specifically, we find steps that occur simultaneouslywith a phenotype transition have distinct sizes comparedwith steps that occur while cells dwell in a particularmorphological phenotype. This makes morphodynamicsa crucial factor in determining the invasive potential ofcancer cells. To our knowledge, this aspect has been sofar largely overlooked in the literature.In the lens of a hidden Markov process, morphody-namics may facilitate cancer invasion because phenotypetransitions allow cancer cells to search for and commit toa more effective migration program [37]. Using a ECMmodel consisting of two mechanically distinct layers, weshow the cells gradually adjust their morphodynamicsas they approach and cross the layer interface (Fig. 6).Therefore morphological phenotype transitions may beessential in cancer cell metastasis by enabling the cells tonavigate non-uniform ECM.In summary, we demonstrate the morphodynamics of3D migrating cancer cells as a powerful tool to inspectthe internal state and microenvironment of the cells. In-vestigated at mesoscale, the morphodynamics imply that3D cancer cell migration is inherently plastic [15]. Theplasticity is controlled by the mode transition matrices,rather than a deterministic decision tree [14]. In order tofurther exploit the information provided by the cell shapefluctuations, future research is needed to decode morpho-dynamics as a rich body language of cells, and to controlmorphodynamics as a route of mechanical programmingof cell phenotype.
MATERIALS AND METHODS
See
SI Appendix
S1-S9 for details of 3D cell culture,microscopy, pharmacological treatments, and data anal-ysis.We thank Prof. Michelle Digman and Prof. StevePress´e for helpful discussions. The funding for this re-search results from a Scialog Program sponsored jointlyby Research Corporation for Science Advancement andthe Gordon and Betty Moore Foundation through agrant to Oregon State University by the Gordon andBetty Moore Foundation. Part of this research was con-ducted at the Northwest Nanotechnology Infrastructure,a National Nanotechnology Coordinated Infrastructuresite at Oregon State University which is supported inpart by the National Science Foundation (grant NNCI-1542101) and Oregon State University. C. Eddy andB. Sun are supported by DOD award W81XWH-20-1-0444 (BC190068). B. Sun is also supported by theNational Institute of General Medical Sciences award1R35GM138179. [1] R. Hooke,
Micrographia (The Royal Society, 1665).[2] P. Mazzarello, Nat. cell biol. , E13 (1999).[3] E. Mayr,
The Growth of the Biological Thought (Belknap,Cambridge, MA, 1982).[4] B. Alberts, A. Johnson, J. Lewis, D. Morgan, M. Raff,K. Roberts, and P. Walter,
Molecular biology of the cell (Garland Science, 2014).[5] R. Singhvi, A. Kumar, G. P. Lopez, G. N. Stephanopou-los, D. I. Wang, G. M. Whitesides, and D. E. Ingber,Science , 696 (1994).[6] M. Diez-Silva, M. Dao, J. Han, C.-T. Lim, and S. Suresh,MRS bulletin / Materials Research Society , 382(2010).[7] A. Serrano-Pozo, M. P. Frosch, E. Masliah, and B. T.Hyman, Cold Spring Harbor Perspectives in Medicine ,a006189 (2011).[8] Z. Yin, A. Sadok, H. Sailem, A. McCarthy, X. Xia, F. Li,M. A. Garcia, L. Evans, A. R. Barr, N. Perrimon, C. J.Marshall, S. T. C. Wong, and C. Bakal, Nat. Cell Biol. , 860 (2013).[9] P.-H. Wu, J. M. Phillip, S. B. Khatau, W.-C. Chen,J. Stirman, S. Rosseel, K. Tschudi, J. Van Patten,M. Wong, S. Gupta, A. S. Baras, J. T. Leek, A. Maitra,and D. Wirtz, Sci. Rep. (2015), 10.1038/srep18437.[10] C. Bakal, J. Aach, G. Church, and N. Perrimon, Science , 1753 (2007).[11] V. K. Lam, T. C. Nguyen, B. M. Chung, G. Nehmetallah,and C. B. Raub, Cytometry Part A , 334 (2017).[12] J. C. Caicedo1, S. Cooper, F. Heigwer, S. Warcha, P. Qiu,C. Molnar, A. S. Vasilevich, J. D. Barry, H. S. Bansal,O. Kraus, M. Wawer11, L. Paavolainen, M. D. Herrmann,M. Rohban, J. Hung, H. Hennig, J. Concannon, I. Smith,P. A. Clemons, S. Singh, P. Rees, P. Horvath, R. G.Linington, and A. E. Carpenter, Nat. Methods , 849 (2017).[13] E. Cai, K. Marchuk, P. Beemiller, C. Beppler, M. G.Rubashkin, V. M. Weaver, A. G´erard, T.-L. Liu, B.-C.Chen, E. Betzig, F. Bartumeus, and M. F. Krummel,Science (2017).[14] C. D. Paul, P. Mistriotis, and K. Konstantopoulos, Na-ture Reviews Cancer , 131 (2017).[15] R. J. Petrie and K. M. Yamada, J. Cell Science , 5917(2012).[16] V. Sanz-Moreno, G. Gadea, J. Ahn, H. Paterson,P. Marra, S. Pinner, E. Sahai, and C. J. Marshall, Cell , 510 (2008).[17] S. Wilkinson, H. F. Paterson, and C. J. Marshall, Nat.Cell Biol. , 255 (2005).[18] V. Sanz-Moreno and C. J. Marshall, Curr. Opin. CellBiol. , 690 (2010).[19] M. Pavlaki and S. Zucker, Cancer Metastasis Rev. ,177 (2003).[20] K. Wolf, I. Mazo, H. Leung, K. Engelke, U. H. von An-drian, E. I. Deryugina, A. Y. Strongin, E. Br¨ocker, andP. Friedl, J. Cell Biol. , 267 (2003).[21] R. J. Petrie, H. M. Harlin, L. T. Korsak, , and K. M.Yamada, J. Cell Biol. , 93 (2017).[22] C. Cortes and V. N. Vapnik, Machine Learning (1995),10.1007/BF00994018.[23] A. Ben-Hur, D. Horn, H. Siegelmann, and V. N. Vapnik,Journal of Machine Learning Research (2001).[24] L. Breiman, Machine learning , 5 (2001).[25] M. Matsubara and M. J. Bissell, Oncotarget ,31602–31622 (2016).[26] A. Daoud, U. Gopal, J. Kaur, and J. S. Isaacs, Oncotar-get , 106807–106819 (2017).[27] T. Yeung, P. C. Georges, L. A. Flanagan, B. Marg,M. Ortiz, M. Funaki, N. Zahir, W. Ming, V. Weaver, andP. A. Janmey, Cell Motil. Cytoskeleton , 24 (2005).[28] L.J.P. van der Maaten, Proceedings of the Twelfth Inter-national Conference on Artificial Intelligence and Statis-tics JMLR W&CP 5 , 384 (2009).[29] G. J. Berman, D. M. Choi, W. Bialek,and J. W. Shaevitz, Journal of The RoyalSociety Interface , 20140672 (2014),https://royalsocietypublishing.org/doi/pdf/10.1098/rsif.2014.0672.[30] K. Yoshlda and T. Soldatl, J. Cell Science , 3833(2006).[31] J. B. Wyckoff, S. E. Pinner, S. Gschmeissner, J. S. Con-deelis, and E. Sahai, Curr. Biol. , 1515 (2006).[32] C. Battle, C. P. Broedersz, N. Fakhri, V. F. Geyer,J. Howard, C. F. Schmidt, and F. C. MacKintosh, Sci-ence , 604 (2016).[33] J.-Y. Tinevez, U. Schulze, G. Salbreux, J. Roensch, J.-F. Joanny, and E. Paluch, Proceedings of the NationalAcademy of Sciences , 18581 (2009).[34] R. J. Petrie, N. Gavara, R. S. Chadwick, and K. M.Yamada, J. Cell Biol. , 439 (2012).[35] M. R. Mejillano, S. Kojima, D. A. Applewhite, F. B.Gertler, T. M. Svitkina, and G. G. Borisy, Cell , 363(2004).[36] A. D. Doyle and K. M. Yamada, Exp. Cell Res. , 60(2016).[37] P. B. Gupta, C. M. Fillmore, G. Jiang, S. D. Shapira,K. Tao, C. Kuperwasser, and E. S. Lander, Cell ,633 (2011). upplementary information for Extracellular Matrix regulates themorphodynamics of 3D migrating cancer cells Christopher Z. Eddy, Helena Raposo, Ryan Wong, Fuxin Li, Bo Sun
October 26, 2020 a r X i v : . [ q - b i o . CB ] O c t
1: Additional experimental details
3D cell culture
GFP-labeled MBA-MB-231 human breast cancer cells are purchased from GenTarget Inc. and are maintainedaccording to the manufacturer’s instructions. To embed the cells in 3D collagen matrices, cells are suspended atvery low density in neutralized collagen solutions. Generally, the suspension is then immediately transferred to glassbottom dishes (ibidi µ -dish) and incubated on either a warming plate set to 25 ◦ C, or in a tissue culture incubator(37 ◦ C, 5% CO2) for 30 minutes in order to solidify the matrix. For fiber alignment, a small magnetic iron shaving(<200 µ m) is first placed onto the dish and then immersed in cell-collagen solution and placed onto the warmingplate. We then drag the particle through the solution by an external magnet along a line for approximately 3 minuteswhile warming, and is then left to solidify (1). The cellularized ECM is then immersed with tissue culture mediumand continuously incubated for 24 hours before imaging. Microscopy and image analysis
Continuous imaging is done with a Leica TCS SPE confocal microscope equipped with a stage-top incubator (ibidi).Images are captured at a rate of 1 frame per 15 minutes. The raw images are gray scale with a resolution of 1024 x1024 pixel . The voxel size has been calibrated to equal 0.538 µ m. A single x-y plane is imaged every 10 µ m in thez-dimension per experiment for up to 24 hours. Using custom Matlab scripts, cell images are maximum-projectedonto a x-y plane and tracked over time (see section S2). The projected images are manually segmented, and screenedto remove cells that are not entirely within the viewing window. S2: Geometric characterization of cell images
Image processing
Following acquisition of fluorescent images, data regarding cell shape and position are obtained by processing andbinarizing the time-lapse images using custom NIH ImageJ and Matlab scripts. First, fluorescence images arebackground subtracted using a rolling ball radius of 50 pixels (26.88 µ m) and then log-transformed in order tomake cell edges highly visible and so that less fluorescent-intense cells are also quantified. A manual threshold isthen applied for each image. After, cells are manually segmented for each z-stack if applicable. Since consecutivez-stacks may have cell overlap, custom Matlab scripts are then used to determine if the same cell is in multiple z-stacks. After, we take a maximum projection (2D) of each cell. Geometrical measurements are then taken on binaryobjects using Matlab’s regionprops function, including area, perimeter, major axis length, minor axis length, solidity,eccentricity, convex area, extent, and equivalent diameter, convex perimeter, fiber length (skeletonized max length),maximum inscribed radius, and maximum bleb radius (maximum secondary circle). Additional measures of formfactor ( AreaP erimeter ), aspect ratio ( MajoraxislengthMinoraxislength ), Convexity (
P erimeterConvexP erimeter ), Curl (
Majoraxislengthfiberlength ), PerimeterCurl (
P erimeterπ (1 − √ − πf ormf actor )), Sphericity ( ∗ MaxInscribedRadiusMajorAxisLength ), Inscribed Area (
MajorAxisLength ∗ πMaxInscribedRadius ),and Bleb Ratio ( MaxBlebRadiusMaxInscribedRadius ) are subsequently calculated, totalling to 21 geometric measures. Collectively,this is shown schematically in figure 1.
Tracking cell position
In order to track the real-space center of the cell, we use maximum inscribed circle (MIC) of the cell image. Com-pared to imaging the nucleus directly, which causes phototoxicity and is prone to poly(nucleus), MIC does notrequire additional probes. To further demonstrate the accuracy of MIC, we compared short videos of dual labeledMDA-MB-231 cells where the GFP channel labels the cytoplasm and RFP channel labels the cell nucleus (SYTO-64,ThermoFisher). We find that MIC agrees very well with direct nucleus staining when determining the cell position,as shown in figure 2. For most of the frames, the deviation is less than 10% of the cell long axis. The root meansquared deviation is approximately 3 microns. 1 ig. 1.
Schematic of image processing and measures taken from binary using custom Matlab and Python scripts(scale-bar = 10 µ m). MEASURE Quantity
Area 490.81 µm Major Axis Length 65.47 µmMinor Axis Length 12.53 µmEccentricity 0.98Convex Area 768.30 µm Equivalent Diameter 25.00 µmSolidity 0.64Extent 0.31Perimeter 167.31 µmConvex Perimeter 146.46 µmFiber Length 73.22 µmMax Inscribed Radius 5.48 µmBleb Length 42.18 µmAspect Ratio 5.23Form Factor 0.02Convexity 1.14Perimeter Curl 10.51 µmCurl 0.89Sphericity 0.17Inscribed Area 2455.99 µmRel. Bleb Length 7.69
Threshold
Measure C oun t s Deviation (µm) Fig. 2.
Positional data acquisition: Position measurements of a cell determined by the nucleus centroid stainedwith SYTO-64 (red) and max-inscribed circle center (green), and histogram of square deviations between the twopositional measures. All scale-bars are 20 µ m. S3: ECM dimension modulates real space and shape space dynamics
The morphodynamics of a cell display a strongly sub-diffusive characterization in shape-space. Shown in Table 1,we quantify the fits of mean square displacement of measures and report the power over a ten hour lag period.2 ∼ x n Measure n (3D) n (2D)Area 0.6261 0.9812Major Axis Length 0.5761 0.8781Minor Axis Length 0.4338 0.4702Eccentricity 0.3638 0.3234Convex Area 0.5738 0.9230Equivalent Diameter 0.5999 0.9246Solidity 0.4493 0.5794Extent 0.4965 0.5208Perimeter 0.5635 0.8353Convex Perimeter 0.6043 0.9636Fiber Length 0.5981 0.8931Maximum-Inscribed Radius 0.5059 0.6059Bleb Length 0.5125 0.6510Aspect Ratio 0.4240 0.6419Form Factor 0.4642 0.5326Convexity 0.2645 0.4477Perimeter Curl 0.4790 0.5681Curl 0.2515 0.2775Sphericity 0.5533 0.5016Inscribed Area 0.5139 0.8617Relative Bleb Length 0.4667 0.6052Real-Space Migration 1.2197 1.3881 Table 1.
Power ( n ) of anomalous diffusion quantified by fitting the mean square displacements of shape and position measures for cellsembedded in a 1.5 mg/mL collagen matrix prepared at room temperature (3D) and cells plated on top of similarly prepared collagenmatrices (2D). S4: Morphological Phenotype Analysis
Details of machine-learning and SVM
In order to classify cells into particular migration mechanisms, we used support-vector machine (SVM) learning (2, 3).As a maximal-margin classifier, SVM was particularly attractive as the overlap between different phenotypes wasunknown. Also important is that cells can display multiple phenotypes at once and thus a soft-margin classifier wasvital. Lastly, given our small dimensional space and small training-sets, SVM was an optimum choice for classificationpurposes.Labeled data were first acquired as described below. Images were binarized and then geometrical data wasobtained on labeled cells for training. We performed parameter grid search for RBF, Linear, and polynomial kernelmodels, with 10-fold cross validation to determine best performance. A grid search determined that an RBF kernelwith γ = 5 . − and C = 5000 yields an average training, validation, and test set accuracy of 93.4%, 86.0%, and94.3%, respectively. However, this model was only nominally improved from the RBF model used with γ = 0 .
01 and C = 1000 (average training, validation, and test set accuracy of 93.0%, 85.6%, and 94.3%, respectively). A linearmodel also recorded strong performance with C = 193 . Details of training, validation, and test sets
Training images were processed as previously described. At least 280 images of each class were acquired to trainmachine learning networks, including the use of an optimized small shearing augmentation (randomly sheared withmagnitude < 0.4) and rotation on images as an over-sampling technique to make our models robust to slight per-turbations and noise, as well as to expand training examples. To evaluate true performance in these networks, anunseen test set was prepared with 35 images. The training and test sets are available on Figshare (4).3 omparisons with Random-Forest classifier
To compare performance to other multi-classification models, we have trained a Random-Forest classifier, tuninghyperparameters of depth, number of estimators, and features (5). The optimized model uses 200 estimators, witha maximum depth of 220 splits, and uses no feature subset selection, yielding a bagged ensemble of regression trees.We perform 10-fold cross validation on 80% of training data and use 20% for validation data, utilizing class weightingto account for the class imbalance during training. We also include our 35 image held out test set for evaluation. Themean accuracy score was 91.3 ± Details of parametric t-SNE embedding algorithm
To visualize a high-dimensional morphological trajectories in three-dimensional space, a t-Distributed StochasticNeighbor Embedding (t-SNE) algorithm is prepared for dimensionality reduction utilizing the geometric character-ization of cells without requiring labels (6). Since the t-SNE algorithm is a non-parametric model, it utilizes allavailable data to separate data clusters well, and cannot classify new data points without retraining and often leadsto different cluster shapes. To embed new data points or entire trajectories in a consistent morphological space,a parametric t-SNE model is prepared. This model uses an Adam optimizer set to minimize the Kullback-Leiblerdivergence loss, and is trained with 800 iterations with a batch size of 256 examples, and the tunable perplexityparameter set to 30.0. The training set consists of 1024 examples, with equal numbers of samples from each pheno-type. The model is made up of three fully connected layers with subsequent ReLU activation after each layer, withan additional final fully connected layer calculating the three t-SNE components for each points. We then use themodel to transform the high-dimensional trajectories of new cells into the three-dimensional space.
S5: Cell Chemical Treatment
To demonstrate the effects of the Rho-ROCK pathway on cellular morphology, we chemically induce activation andinhibition of pathway proteins. Y27632, a specific inhibitor of Rho-associated protein kinases as well as ROCK-IIactivity (7–9), was purchased from Sigma-Aldrich and diluted to a working concentration at 3 µ g/mL (0.1% v/vDMSO) in serum-free growth medium. Similarly, Rho activator II (CN03), known to robustly increase the level ofGTP bound RhoA (10, 11), was purchased from Cytoskeleton and diluted to a working concentration at 2 µ g/mL(0.1% v/v DMSO) in serum-free growth medium. 1.5 mg/mL gels were prepared as discussed in the main text withcells at a low concentration and incubated in serum-free growth medium for 6 hours. Cells were then imaged withconfocal microscopy for up to 4 hours. Then, serum-free growth medium is removed from the dish and replaced withfresh serum-free growth medium containing the prepared chemical. We then proceed to immediately image for upto another 20 hours.Figure 4 shows additional trajectories of cells that responded to chemical treatment of CN03. Generally, cellstreated with CN03 exhibit morphologies characteristic of cell contraction. Notably, cells do not produce protrusionsfollowing treatment and typically retract protrusions post-treatment. Visually, this is seen as a notable slide towardthe blebbing migration mode (green) or toward lamellipodial cell-spreading (blue).Conversely, figure 3 shows additional trajectories of cells that responded to chemical treatment of Y27632. Cellstreated with Y27632 show characteristic production of protrusions and/or sustained pre-existing protrusions. In thefigure, this is shown by most trajectories moving toward filopodial migration (magenta) or toward actin-enrichedleading edge migration (yellow). 4 ig. 3. Additional t-SNE time projections (solid) of ROCK-inhibiting Y27632 drug treated cells.
Fig. 4.
Additional t-SNE time projections (solid) of Rho-activating CN03 treated cells.
S6: ECM Characterization
Gel Autocorrelation
In order to quantify the density fluctuations of collagen ECM, we compute the autocorrelation functions of confocalreflection images of collagen gels. Reflection images are first background subtracted using a rolling ball with radius5f 50 pixels (26.88 µ m). Images were then log-transformed to make fibers highly quantifiable. Images are then meansubtracted, and the autocorrelation is calculated. Following, the autocorrelation is normalized and then smoothedusing interpolation. The results are shown in Figure 5. The spatial uniformity indicates an appropriate distributionof randomly oriented gel fibers. Additionally, these show that the decay in the autocorrelation is slower for higherdensity gels, an indication of the fiber quantity, and much faster for higher temperature gel, caused by the muchshorter fiber lengths and smaller pores (12). The anisotropy of decay in the final graph indicates the direction ofalignment. M agn i t ude M agn i t ude
002 2 1 M agn i t ude
02 022022 Distance (µm) D i s t an c e ( µ m ) Distance (µm) D i s t an c e ( µ m ) Distance (µm) D i s t an c e ( µ m ) M agn i t ude
02 2022 Distance (µm) D i s t an c e ( µ m ) Fig. 5.
Autocorrelations shown for 25 ◦ C 1.5 mg/mL Randomly oriented gel (upper left), 25 ◦ C 3.0 mg/mL Randomlyoriented gel (upper right), 37 ◦ C 1.5 mg/mL Randomly oriented gel (lower left), and 25 ◦ C 1.5 mg/mL Aligned orientedgel (lower right).
Gel Coherency
In order to determine the degree of local and global alignment of collagen fibers, we take the pre-processed confocalreflection images and use OrientationJ with 9-14 circular ROIS per image (packed without overlap), with ROI sizesranging from 145 to 450 pixel diameter. At least 3 images (one per experiment) were used to quantify coherency.Empirically, larger ROIs were reliable to quantify global alignment, while smaller ROIs were better used for localalignment measures. Smaller ROIs than 145 pixels were not used, as it was noticed the coherency measured byOrientationJ can be highly biased for large fibers such as in 25 ◦ C gels (145 pixels [ 78 µ m] is approximately twicethe average fiber length [ 41 µ m] in a 25 ◦ C gel). The results are shown in figure 6 and follows previous literature(13). C ohe r en cy ** RandomAligned C ohe r en cy Fig. 6.
Coherency measurements taken using custom Matlab scripts utilizing OrientationJ plug-in for ImageJ (NIH).Measurements are for 1.5 mg/mL collagen gels prepared at room temperature, where we have utilized the alignmentprotocol to compare against typical randomly-oriented collagen gels.6 heology
Strain sweep rheology measurements were performed on the varying density gels with a Discovery Hybrid Rheometer-3 (TA instruments) at a 1 Hz frequency in a parallel plate geometry. A standard peltier plate (TA instruments)allowed gels to be formed at 37 ◦ C. Young’s modulus is shown in figure 7. As shown, storage moduli in the linearregime for gels with 1.5, 3.0, and 4.5 mg/mL collagen were about 70, 180, and 250 Pa, respectively. -1 Strain (%) G ' ( P a ) Fig. 7.
Strain-sweep rheology measurements of collagen fiber networks gelled at 37 ◦ C using DHR-3 Rheometer (TAInstruments) in a parallel plate geometry at a 1 Hz frequency.
S7: Additional details of morphological phenotype dynamics
Dwell time definition
The dwell time of state i is determined by using D i → i = 11 − r i → i . where r i → i is the transition rate from state i back to itself, and D i → i is the corresponding dwell time. We found thisdefinition of dwell time maximizes the usefulness of data and is more robust to experimental shortfalls that affectthe naive method of simply counting the number of frames a cell remains in the same phenotype for. Details of transition rate calculations
Assuming the probability distribution of cell morphological phenotypes follow a Boltzmann distribution, then theprobability, P i , given by P i = 11 + P Nj = i r i → j r j → i where r i → j is the transition rate from state i to state j , and conversely r j → i is the transition rate from state j tostate i . The transition matrix is thus calculated as probability per unit time, and hence each row in the transitionmatrix will sum to 1 ( r i → i −
1) + N X j r i → j = 1 . The transition rate r i → j is simply calculated by counting the number of transitions from state i to state j , and thendividing by all observed transitions from state i . Using this method, we find that the transition rates are stable,7egardless of the length of trajectories in computational experiments.Following SVM classification, phenotype dynamics could be properly drawn from data. Importantly, where max-imum decisions by the SVM classifier do not exceed 60%, the classification is thus determined to an intermediatebetween two states. The intermediate state can be a chimera of two states, with most occurrences being as a celltransitions between morphological phenotypes. Because this state is not considered to be unique, we calculate thetransition rate from state i that passes through N intermediate states prior to state j as 1 /N frames . We find bysimulation that this method can most-accurately recover transition rates in comparison to methods using soft-maxor ignoring intermediate states. Probability flux calculations
To investigate broken detailed balance in morphological phenotype space, we report probability flux calculationsshown in figure 8. The probability flux from state i to state j is given by∆ φ i → j = N i → j P k P l N k → l where N i → j is the number of transitions from state i to state j , and both summations in the denominator are overall available states. Therefore, the net probability flux between states i and j is given by ∆ φ i → j − ∆ φ j → i . We findthat the net flux between states are all less than 0.003 probability difference per hour for cells in any ECM conditionwe tested. To reveal if any small difference in probability flux may be significant, we also report the probability fluxpercent difference between states i and j is given by∆ φ i → j − ∆ φ j → i ∆ φ i → j + ∆ φ j → i We find that the maximum net probability flux percent difference for cells in any ECM condition evaluated is lessthan 9%.
A BC D prob / hr AE FP BB LA
E FG H % difference 9.00
AE FP BB LA
Fig. 8. (A-D) Net probability flux calculations under varying ECM conditions. (E-H) Percent difference of netprobability flux under varying ECM conditions. (A, E) collagen ECM prepared at room temperature (RT) and [ col ]= 1.5 mg/mL. (B, F) collagen ECM prepared at RT and [ col ] = 3.0 mg/mL. (C, G) collagen ECM prepared at 37 ◦ C and [ col ] = 1.5 mg/mL. (D, H) collagen ECM prepared at RT and [ col ] = 1.5 mg/mL with flow-aligned collagenfibers. 8
8: Motility analysis of morphological phenotypes
In order to determine the migrational persistence of migration modes, we report first the velocity autocorrelationobserved for cells in collagen ECM prepared at room temperature (RT) and [ col ] = 1.5 mg/mL. We find that theautocorrelation quickly decays to zero in a single time step (15 minutes) shown in Figure 9(A). We additionallychecked the cos ( θ ) distribution between consecutive steps and found steps that the concurrent steps seem to betaken randomly in direction, as given by the U-shaped distribution of Figure 9(B).
0 1 2 3 4 5
Time (hr) -0.200.20.40.60.81 VA C F ( t ) A B
Fig. 9. (A) Velocity autocorrelation and (B) binned directionality given by cos ( θ ) between consecutive velocityvectors for cell migration in collagen ECM prepared at room temperature and [ col ] = 1.5 mg/mL. (Bin width = 0.1)A soft-max classification strategy was used to determine step-sizes dominated by a single migration mode andto avoid confounding step sizes from classifications near the decision function boundary. To evaluate the motilityof various migration mode transitions, step sizes were separated into components parallel and perpendicular to thedirection of the previous step. Both directions had step-size distributions that were approximately Gaussian. Figure10(A) shows the binned step size distribution in the persistent direction fit with a Gaussian distribution. We findthat these fits yield mean step sizes close to zero and small variances, which are not concurrent with the mean-squaredisplacements (MSD) of RT data shown in Figure 5B. Rather, we find more suitable fits are obtained by fitting themagnitude of the step sizes, which follow a log-normal distribution as shown in figure 10(B). These fits result in stepdistributions that yield similar MSD values to RT data. The parameters used to generate the log-normal fits arecalculated from fitting the empirical cumulative distribution of step size magnitudes with a log-normal CDF, shownin figure 11. 9 B Fig. 10. (A) Gaussian fitting (red) of binned step sizes of in the persistent direction (blue) [bin width = 1 µ m] forvarious migration mode transitions. (B) Log-normal fitting (red) of binned step size magnitudes (blue) using theparameters shown in Figure 5A to generate the log-normal PDF [bin-width = 0.5 µ m] for various migration modetransition. Step Size ( m) C u m . P r ob AM-AM
Step Size ( m) C u m . P r ob AM-ME
Step Size ( m) C u m . P r ob ME-AM
Step Size ( m) C u m . P r ob ME-ME
Fig. 11.
Empirical cumulative distribution (blue) and log-normal cumulative distribution function fit (red) of stepsizes for transitions between amoeboid (AM) and mesenchymal (ME) modes.
S9: Interface calculations
Experiment Details
Interface experiment was done in triplet. Briefly, the outer gel was first made by gelling cells in collagen solution (1.5mg/mL neutralized) at 25 ◦ C for 20 minutes on the DIGME stage (14), warmed at 25 ◦ C. After, the needle was gentlyremoved and a new ice-cold collagen solution (1.5 mg/mL neutralized) containing cells was then dripped into the hole,gently swirled, and then placed into the incubator at 37 ◦ C for 15 minutes. After, 3 mL of growth medium was addedon top and sat for 24 hours before imaging. Imaging was taken near the interface, imaging bright-field, fluorescence(green), and back-reflection confocal images. Using the back-reflection images, the interface was manually traced outthrough the z-stack. The distance to the interface from cell centroids in corresponding z-stacks were then measuredaway from the closest marked interface point using Matlab bwdist. The frequency vs distance is shown in Figure 12,indicating a large number of cells were images close to the interface.10 istance (µm) C oun t s Fig. 12.
Frequency of cell locations away from interface (at 0). Negative is in 37 ◦ C gel, Positive indicates 25 ◦ C gel).
Details of interface analysis
Migration mode transitions were first determined by using continuous trajectories accounting properly for interme-diate state classification, and then distances away from the interface were determined by the initial state location.A 1-D Gaussian kernel was used to acquire continuous local spatial probability densities of transitions (per hour),centered every 5.376 µ m (10x the distance-to-pixel ratio) with a standard deviation of 26.88 µ m. This yields theprobability density (per hour) of observing a transition P ( i → j, x ) at a location x, as used in the main text. Tomitigate the bias from non-uniform cell density, we then divide the prior probability by the spatial probability ofobserving a cell within the spatial window centered at x, M ( x ). M ( x ) is calculated using the same Gaussian kernel. [1] Jihan Kim, Yu Zheng, Amani A. Alobaidi, Hanqing Nan, Jianxiang Tian, Yang Jiao, and Bo Sun. Geometric dependence of 3d collectivecancer invasion. Biophysical Journal , 118(5):1177 – 1182, 2020. ISSN 0006-3495. doi: https://doi.org/10.1016/j.bpj.2020.01.008. URL .[2] Corinna Cortes and Vladimir N. Vapnik. Support-vector networks.
Machine Learning , 20(3), 1995. doi: 10.1007/BF00994018.[3] Asa Ben-Hur, David Horn, Hava T. Siegelmann, and Vladimir Vapnik. Support vector clustering.
J. Mach. Learn. Res. , 2:125–137, March 2002.ISSN 1532-4435.[4] Christopher Eddy. MDA-MB-231 Morphology Training and Test Sets.
FigShare , 10 2020. doi: 10.6084/10.6084/m9.figshare.13125086.[5] Leo Breiman. Random forests.
Machine learning , 45(1):5–32, 2001.[6] Gordon J. Berman, Daniel M. Choi, William Bialek, and Joshua W. Shaevitz. Mapping the stereotyped behaviour of freelymoving fruit flies.
Journal of The Royal Society Interface , 11(99):20140672, 2014. doi: 10.1098/rsif.2014.0672. URL https://royalsocietypublishing.org/doi/abs/10.1098/rsif.2014.0672 .[7] Renaud Poincloux, Olivier Collin, Floria Lizárraga, Maryse Romao, Marcel Debray, Matthieu Piel, and Philippe Chavrier. Contractility ofthe cell rear drives invasion of breast tumor cells in 3d matrigel.
Proceedings of the National Academy of Sciences , 108(5), 2011. doi:10.1073/pnas.1010396108.[8] Simon Gordonov, Mun Kyung Hwang, Alan Wells, Frank B. Gertler, Douglas A. Lauffenburger, and Mark Bathe. Time series modeling of live-cellshape dynamics for image-based phenotypic profiling.
Integrative Biology , 8(1), 2015. doi: 10.1039/c5ib00283d.[9] Masahiro Matsubara and Mina J Bissell. Inhibitors of rho kinase (rock) signaling revert the malignant phenotype of breast cancer cells in 3d context.
Oncotarget , 7(22), 2016. doi: 10.18632/oncotarget.9395.[10] Meng-Horng Lee, Pei-Hsun Wu, Daniele Gilkes, Ivie Aifuwa, and Denis Wirtz. Normal mammary epithelial cells promote carcinoma basementmembrane invasion by inducing microtubule-rich protrusions.
Oncotarget , 6(32):32634–32645, 2015. ISSN 1949-2553. doi: 10.18632/oncotar-get.4728.[11] Abdelkader Daoud, Udhayakumar Gopal, Jasmine Kaur, and Jennifer S Isaacs. Molecular and functional crosstalk between extracellular hsp90and ephrin a1 signaling.
Oncotarget , 8(63), 2017. doi: 10.18632/oncotarget.22370.[12] Christopher Allen Rucksack Jones, Long Liang, Daniel Lin, Yang Jiao, and Bo Sun. The spatial-temporal characteristics of type i collagen-basedextracellular matrix.
Soft Matter , 10(44):8855–8863, November 2014. ISSN 1744-683X. doi: 10.1039/c4sm01772b.[13] T. D. Clemons, M. Bradshaw, P. Toshniwal, N. Chaudhari, A. W. Stevenson, J. Lynch, M. W. Fear, F. M. Wood, and K. Swaminathan Iyer. Coherencyimage analysis to quantify collagen architecture: implications in scar assessment.
RSC Adv. , 8:9661–9669, 2018. doi: 10.1039/C7RA12693J.URL http://dx.doi.org/10.1039/C7RA12693J .[14] Amani A. Alobaidi and Bo Sun. Probing three-dimensional collective cancer invasion with digme.
Cancer Convergence , 1(1), 2017. doi:10.1186/s41236-017-0004-9., 1(1), 2017. doi:10.1186/s41236-017-0004-9.