Lévy-like movement patterns of metastatic cancer cells revealed in microfabricated systems and implicated in vivo

Metastatic cancer cells differ from their non-metastatic counterparts not only in terms of molecular composition and genetics, but also by the very strategy they employ for locomotion. Here, we analyzed large-scale statistics for cells migrating on linear microtracks to show that metastatic cancer cells follow a qualitatively different movement strategy than their non-invasive counterparts. The trajectories of metastatic cells display clusters of small steps that are interspersed with long “flights”. Such movements are characterized by heavy-tailed, truncated power law distributions of persistence times and are consistent with the Lévy walks that are also often employed by animal predators searching for scarce prey or food sources. In contrast, non-metastatic cancerous cells perform simple diffusive movements. These findings are supported by preliminary experiments with cancer cells migrating away from primary tumors in vivo. The use of chemical inhibitors targeting actin-binding proteins allows for “reprogramming” the Lévy walks into either diffusive or ballistic movements.

T he motility of mammalian cells has been studied for decades 1,2 , and trajectories of cell movements have been quantified in various ways. Early models of cell motility were founded on the classic Langevin equation and described the movements of adherent cells [3][4][5] (for description of smaller, faster, and weakly-adherent immune cells, see ref. 6,7 ) as an Ornstein-Uhlenbeck (OU) process, 8 such that the cell's mean square displacement, < x(t) 2 > , is expressed as 2nD (t -P (1exp (-t/P)), where n is the dimensionality of the system, D is the diffusion coefficient, and P is the so-called persistence time. This model predicts Gaussian distribution of velocities that are exponentially correlated in time, leading to directional persistence on short time scales (t < P) 9 . On long time scales (t » P), the model is reduced to a random walk and predicts uniform distribution of turning angles. This formalism has been successfully applied to describe the motions-coined persistent random walks (PRWs)of fibroblasts, lung epithelial cells, and microvessel endothelial cells [3][4][5] . In addition, due to apparent directional persistence, animal as well as cellular movements have frequently been described as "correlated random walks" (CRWs) 7,10,11 . In CRWs 1 , the step sizes/times are drawn from the Gaussian or other exponentially decaying distribution, and the direction of the preceding step influences the direction of the next step. Overall, both PRW and CRW models predict that cell movements are ballistic (i.e., persistent in direction or < x 2 >~t α with α = 2) at short times and diffusive (α = 1) at long time scales. Lévy walksrecently detected in T cells 6 -are different because they are superdiffusive 12 (i.e., < x 2 >~t α with 1 < α < 2) at all times and are composed of sequences of many short steps interspersed with longer "flights". This pattern is conserved across all scales, in effect giving rise to fractal patterns with no characteristic scale 13 . Mathematically, Lévy walks 14,15 are characterized by non-Gaussian, heavy-tailed, power law distribution of persistence times/ step sizes, such that P(t)~t -μ , where t is persistence time/step size or time/distance it takes to move one step between the turns and μ is power law (Lévy) exponent with 1 < μ < 3. In the absence of a characteristic scale, the overall length of a Lévy walk is determined by the longest step and the step-length variance grows over time, though it remains finite even when unbounded by biological and environmental considerations. Most biological systems are bounded/limited (e.g., cell trajectories are limited by cell cycle and environmental conditions), resulting in truncation of the power law tail 14,16-18, which introduces characteristic scale to the movement pattern. However, variability around the characteristic scale is very large and self-similar, which is in sharp contrast to other finite-scale movement patterns. Lévy-like movement patterns have been observed in movements of a number of multicellular animals [17][18][19][20][21][22][23][24][25] and humans 26,27 , found in trace fossil trails 17 , and recently also in T cells searching for parasite-infected cells 6 , swarming bacteria 28,29 , and even molecular motors within cells 16 . The observation of Lévy walks has been attributed to the execution of an optimal search strategy for sparsely and randomly distributed resources/target sites [30][31][32][33] , though this interpretation has not been generally accepted 13,34 .
With the general applicability of each of these models still being debated, 9,35 one outstanding and important question concerns the movement patterns of non-metastatic versus metastatic cancer cells. Although the latter are known to have higher migration velocity and, in some cases, increased directional persistence [36][37][38][39] , it remains unclear whether increased metastatic potential is reflected by a qualitative change in the overall strategy of cell locomotion.
In this work, we address this question with a material system of micropatterned lines 10,40 on which the cells perform onedimensional motions which (i) have been shown 41 to mimic the motility of cells migrating in 3D better than commonly used planar substrates and (ii) are observed in metastasizing tumor cells in vivo, attaching preferentially to and moving along linear fibers (e.g., collagen fibers 36 ) or along preexisting linear perimuscular or perineural "microchannels" 42,43 . Importantly, the "microtracks" we use enable unambiguous determination of persistence times-as the times that cells move "to the left" or "to the right" before reverting their direction of motion 44 -and collection of large numbers of data points such that even the lowprobability events are captured. Statistical analysis of such data then reveals that while non-metastatic cancerous cells are simple diffusive movers, the metastatic ones not only move in a superdiffusive fashion (which has been shown before 45,46 but not in the context of metastatic cells), but also display movement patterns consistent with Lévy walks 1,47 . Significantly, we also detect Lévy walks of metastatic cells migrating in preliminary in vivo studies, in which we used state-of-the-art intravital multiphoton microscopy to resolve trajectories of individual cancer cells migrating away from animal-implanted tumors. While the generality of the observed trends certainly merits additional studies (especially in vivo), the results we describe pose some intriguing questions as to why the invasive cancerous cells have developed a movement strategy that is frequently observed in multicellular animal predators and hunters [17][18][19][20][21][22][23][24][25] (as well as killer/effector T cells 6 searching for rare target cells) and, as mentioned above, is thought to correspond to an optimal search strategy for sparsely and randomly distributed resources/target sites [30][31][32][33] (but see also ,34 ). In this context, toward the end of the paper, we describe experiments in which RNA interference and chemical inhibition of actin-binding proteins can change the Lévy walking phenotype of metastatic "predators" into either unidirectional, ballistic motions, or into diffusive migrations characterizing benign or non-invasive cancerous cells.

Results
Lévy walks of metastatic cancer cells revealed on linear microtracks. Most of the in vitro experiments were performed on linear microtracks etched in gold-on-glass substrates using the socalled Wet Etching technique 40,48,49 (Fig. 1a; for all fabrication details, see Supplementary Methods/Supplementary Note 1). The gold regions were protected with self-assembled monolayers, SAMs 50,51 , of oligo (ethylene glycol)-terminated alkyl thiols (HS (CH 2 ) 11 (OCH 2 CH 2 ) 6 OH; EG 6 (ProChimia Surfaces, Gdansk, Poland: www.prochimia.com) known to prevent cell adhesion. The unprotected glass lines were covered with either Laminin (Sigma-Aldrich, cat. # L2020) or Laminin 5 (LN 5; extracted from 804G cells in crude form, as previously described 52 )-these extracellular matrix substrates were chosen because of their motility-promoting characteristics (vs. more adhesive fibronectin) and because they are routinely used as physiologically relevant substrates of choice for the respective cell lines.
When the cells were applied (at plating density of~10,000 cells/cm 2 ) onto microstructured substrates presenting arrays of 20-μm-wide linear tracks, they localized exclusively onto these tracks, spread, and, to a good approximation, displayed onedimensional motions (Fig. 1b). We compared and contrasted motions of six types of cells from three cancers ( Fig. 2; Supplementary Figure 1): non-metastatic PC-3 and metastatic PC-3M 53 prostate cancer cells; non-metastatic MCF-7 and metastatic MDA-MB-231 38 breast cancer cells; and nonmetastatic B16-F0 and metastatic B16-F1 54 mouse melanoma cells. Regarding the cell line choices, we note that for B16 and PC lines, cells are termed metastatic versus non-metastatic based on, respectively, their low and high metastatic potentials 53,54 . For breast cancer lines, the MCF-7 cell line retains several characteristics of differentiated mammary epithelium and represents a poorly invasive luminal subtype of breast cancer, whereas the MDA-MB-231 line represents a highly invasive basal subtype of breast cancer 55 .
Using an automated image acquisition and analysis protocol developed in-house, we were able to monitor cell motions for up to 16 h and collect robust statistics with n = 17-69 different cells and~5120-20,800 time points per every cell type studied (see Table 1 legend). Only tracks housing one cell (~60% of all tracks at the plating density used) were analyzed to eliminate any potential artifacts due to cell-cell collisions (see Supplementary Note 1 for comment on cell collisions). The typical trajectories of cells on the microtracks shown in Fig. 1c, Fig. 2a To quantify these cell trajectories and the dynamics of cell motions, we first plotted cells' mean square displacements, , versus time, t, where ::: h i denotes a combined average over all starting times t 0 and cell paths (see refs. 12,45 ). When plotted on a log-log scale, these dependencies give straight lines corresponding to x 2 h i ¼ t α scaling, with the slope of the lines being the exponent α. From statistical physics, it is well known 12 that α = 1 indicates diffusive motion, values 1 < α <2 correspond to superdiffusion, and α = 2 is characteristic of ballistic motion (e.g., unidirectional motion without turns). As evidenced by the plots in Fig. 2c, all non-metastatic cells we studied move diffusively (α~1) while their highly metastatic variants move in a superdiffusive fashion (α~1.5-1.6). These characteristics are conserved over the entire time domain and not only short term (as CRW or PRWs).
Long-term superdiffusive behavior is sometimes found in systems performing Lévy walks, 47 which are a sub-class of random walks and are characterized by periods of small steps interspersed with long but infrequent unidirectional excursions. Typical cell trajectories such as those shown in Fig. 1c and Fig. 2a [13][14][15] suggest that the metastatic MDA-MB-231, PC-3M, and B16-F1 cells might indeed be performing Lévy walks. This hypothesis can be mathematically verified by the analysis of cells' cumulative frequency distributions, CFDs, of persistence times (i.e., the persistence time is defined as the time that cell moves persistently in a given direction 12 ). When such distributions are plotted on a log-log scale, it is evident that the CFDs corresponding to metastatic cells are more heavy-tailed than CFDs corresponding to non-metastatic ones ( Fig. 2 d-f). However, extreme care must be taken when assigning such trends to a specific functional form, especially given recent findings 56,57 that the use of inaccurate statistical methods has led to an incorrect assignment of the search/movement patterns as Lévy flights in a number of studies. Accordingly, we followed the rigorous procedure developed by Edwards et al. 57,58 to test power law distributions using the likelihood and Akaike weights (which measure an appropriateness of a given fit 57,58 ) and tested multiple alternative competing models (see also Supplementary Notes 1 and 2, and Supplementary Tables 1, 2) as detailed by Clauset et al 56 . In addition to power law and exponential models, we also considered heavy-tailed log-normal (observed in motions of T cells within lymph nodes 7 ), heavy-tailed stretched exponential 56 , as well as a truncated power law 59 (in which power law tail of Lévy walk is truncated).
In particular, the truncated power law was considered because data collection time was limited for both short and long times. The short time limitation at~3-5 min was due to relatively large cell sizes and small speeds (respectively,~20-40 μm and~0.25-2 μm/min for metastatic cells vs.~7 μm and~6-12 μm/min for T cells studied in ref. 6 ), making it difficult to (1) distinguish their translocation from just the membrane dynamics, and (2) limit the photo-damage imposed to cells by laser/light exposure. The long time limitation was due to the cells dividing (e.g.,~30 h for MDA-MB-231). It is therefore impossible to obtain single-cell data spanning the desirable two decades of persistence times and enhancing statistics of lowprobability, long cellular excursions in the power law's tail 60 . Instead, we (akin to others 9,35,46   The results summarized in Fig. 2d-f and Table 1 demonstrate that only the metastatic cells (PC-3M, MDA-MB-231, B16-F1) exhibit log-log distributions of persistence times that can be assigned as Lévy walks. While the cumulative frequency distribution, CFD, plots might "appear" to fit a simple power law (PðtÞ ¼ ðμ À 1Þa μÀ1 t Àμ ; t ! a, where ðμ À 1Þa μÀ1 is the normalization constant and parameter a determines the time at which a power law fit is appropriate; a values are listed in Table 1), detailed maximum likelihood and Akaike-weight analyses evidence that a more appropriate fit is a truncated power law, P t where b is an upper truncation parameter corresponding to the maximal point in each dataset. Importantly, the values of exponents µ in these fits are below 3, corresponding to Lévy walks. Moreover, such analyses confirm that the truncated power law fit is superior to all alternative models considered ( Fig. 2d- Table 1 ARTICLE NATURE COMMUNICATIONS | DOI: 10.1038/s41467-018-06563-w and 2 and Supplementary Methods for all additional details on model fitting and comparisons). All data sets for metastatic cells pass the goodness-of-fit test performed according to the method detailed in ref. 56 Figure 7). In this context, it is important to note that truncated power law (with µ~1-3) has become increasingly favored (over pure power law) as a more relevant model for describing bounded/limited biological systems and-similar to a pure power law-is considered Lévy walk 14,[16][17][18] . We also note that (i) "walks" rather than "flights" 15 terminology is appropriate since, as we verified, our migrating cells exhibit strong correlation between persistence times and persistence lengths (i.e., step sizes); and (ii) observing Lévy-like patterns in a simple 1D system in the absence of chemoattractant gradients suggests that they may be an inherent characteristic of metastatic cells.

(Supplementary Note 3 and Supplementary
Regarding non-metastatic cells, their CFDs also fit best to a truncated power law, but the exponents μ in these fits are greater than 3 (Fig. 2d-f and Table 1), which is in line with the purely diffusive 34,47 nature of these cells' motions.
Lévy walks of metastatic cancer cells in live tumors. While the microtracks offer a convenient platform for collecting large numbers of motility data, these ex vivo results are not necessarily indicative of in vivo motility patterns, which may be influenced by local chemotactic gradients and/or cell interactions with the tissue microenvironment. Extension to in vivo, however, is technically challenging 36,61 and further complicated in our case where resolving individual cells over long (hours) times is necessary for the collection of statistics ample enough to construct reliable probability distributions. Also, it should be remembered that in vivo analyses cannot offer complete correspondence with all the cell types studied on microtracks and are limited to tumors (i) lying within the penetration depth of the existing high-resolution microscopy modalities (max. 600 μm in mouse dermis), and (ii) for which appropriate animal protocols have been validated and approved. In our study, the additional requirement to compare metastatic versus non-metastatic cells of the same origin limited the available choices to a pair of non-metastatic B16-F0 and metastatic B16-F10 mouse melanoma cells. The trajectories of these cells invading the mouse dermis were recorded using a dorsal skin-fold chamber model 43 . Using combined near-infrared and infrared multiphoton microscopy 43 , second-harmonic generation (SHG) allowed for the reconstruction of fibrillar collagen and myofibers, while fluorescence enabled visualization of blood vessels and tracking of moving cell nuclei (Fig. 3).
On one hand, both B16-F0 and B16-F10 cells predominantly formed invasion strands consisting of up to 100 individual cells which jointly infiltrated the dermis by following linear tissue interfaces provided by parallel perimuscular and perineural microtracks, blood vessels, or collagen bundles (Fig. 3). On the other hand, the path organization and kinetics of cell motions  Table 2). The path organizations were found to be strongly zonedependent. At the tumor margin (zone 1), where cells are "crowded" and frequent cell-cell collisions occur, the motions of both B16-F0 and B16-F10 cells are "jiggly" and purely diffusive, as evidenced by the "diffusion exponent" α~1 (Fig. 4c). Beyond the tumor border, within zone 2, the cells experience more frequent collisions from one side (zone 1) and are effectively "pushed away" from the tumor exhibiting slightly differing degrees of superdiffusive motions (α~1.38 for B16-F0 and α~1.63 for B16-F10). The differences between the inherent characteristics of cell motions-i.e., not dominated by crowding effects/cell-cell collisions-manifest themselves fully in zones 3 and 4, with nonmetastatic cells reverting to diffusive motions (α~1) while metastatic cells remaining evidently superdiffusive, with α1 .47-1.87. These characteristics were conserved over the entire time domain (i.e., not only short term as in correlated/persistent random walks) and were statistically significant as evidenced by the 95% confidence intervals for α listed in Table 2.
Most importantly, the cumulative frequency distributions, CFDs, of persistence times shown in Fig. 4d reveal that nonmetastatic B16-F0 cells move diffusively (truncated power law with μ > 3), whereas the metastatic B16-F10 cells have a CFD characteristic of a Lévy walk (truncated power law with 2 < μ < 3). For the B16-F0 cells, the Akaike weights quantifying the appropriateness of the fits show that even though truncated power law fit is preferred over pure power law or exponential fits, μ > 3 indicates diffusive motion for all zones. For B16-F10 cells, truncated power law fit is also appropriate in all zones with Akaike weights close to 1, but the value of exponent 2 < μ < 3 for zones 2-4 indicates Lévy-like motion ( Table 2). When analyzing all data for B16-F10 from zones 2-4 together, µ~2.36 ± 0.13 with Akaike weights for truncated power law close to 1. Overall, these in vivo studies resemble the results obtained for cell's on microtracks both in terms of the diffusive-vs-Lévy dichotomy and the linearity of cell trajectories (here, along tissue microchannels rather than microtracks). We note that additional analysis of cells migrating through 3D collagen gels (to probe the effects of 1D vs. 3D microenvironment) are also consistent with diffusive-vs-Lévy motions for normal versus cancer cells (see Supplementary Note 4 and Supplementary Figures 8, 9).
"Reprogramming" Lévy walks into other types of motions. Assuming that the Lévy-type walks might be an untoward characteristic of "predatory" metastatic cells, we next focused on the question whether these walks could somehow be altered-in particular-could they be reverted into diffusive motions characterizing the non-invasive cancer cells?
To this end, we inhibited-either chemically or by using specific siRNAs-selected proteins known to be involved in polymerization of actin filaments and their organization into higher-order structures and known to drive extension of cellular protrusions 62 . Results summarized in Fig. 5 and Table 3 (see also Supplementary Movies 9-12 and 16-18) for the metastatic MDA-MB-231 cells moving on the microtracks reveal that inhibiting Rac1 (with chemical inhibitor NSC23766) or knocking down of Cofilin-1 (involved in actin filament depolymerization), Profilin-1 (regulating actin polymerization), or Dia-1 (involved in polymerization of unbranched actin filaments) increased μ exponents slightly, but not as much as to fully eliminate Lévy-like motions, (Table 3, see also Supplementary Figures 10-13). In contrast, inhibition of the Arp2/3 (but not of upstream GTPase Rac1, Fig. 5b and Supplementary Movie 11)-involved in the polymerization of actin into dendritic/branched networks at cell's "front"-by a small-molecule inhibitor CK666 not only slowed the metastatic cells down, from 0.99 μm/min to 0.57 μm/min, but  Table 2 changed the nature of their motions from Lévy-like to diffusive (Fig. 5d, Supplementary Figure 14, Supplementary Movies 9, 17, and Supplementary Table 4 for the summary of all speed data). When Myosin II-typically involved in forming contractile actin networks and bundles and causing concurrent cell rear retraction 62 -was inhibited by blebbistatin at intermediate inhibitor concentrations (see Fig. 5 legend), the cell speed changed only slightly from 0.99 μm/min to 1.13 μm/min, but caused cells to move ballistically (i.e., steadily in one direction; note that in the context of this work, "ballistic" corresponds to motions for which persistence times are longer than total observation time and for which α~2; Fig. 5c and Supplementary Movie 10). In addition, pronounced changes in the motility patterns could also be achieved by inhibition of several proteins simultaneously. A case in point here is the reversal of Lévy to diffusive walks by the simultaneous inhibition of Myosin II with 10 μM of blebbistatin and of Rac-1 with 100 μM of NSC23766 ( Supplementary  Figures 11 and 13, Supplementary Movies 12 and 18)-what is interesting about this result is that the inhibition of each of the proteins separately (cf. Fig. 5b, c) does not lead to diffusive walking, implying a synergistic/cooperative action of the two proteins in controlling the Lévy-walk movement pattern.
The role of protrusion-retraction synchronization in Lévy walks. Our previous morphodynamic profiling analyses have shown that in metastatic cells, the protrusions and retractions are highly "synchronized" both in space and in time; in contrast, protrusions and retractions formed by non-metastatic cells are not "synchronized" 63 . In order to test if such protrusionretraction synchronization can give rise to truncated Lévy walk/ power law step-size distributions, we developed a simple model described in detail in Supplementary Note 5, Supplementary  Figures 15 and 16. Briefly, we considered a scenario in which the probabilities of taking left-right steps depend on prior history (non-Markovian process). We showed that if consecutive steps slightly favor moves in the "same" directione.g., as observed in the microscale dynamics of cell membrane where persistence of membrane protrusion typically depends on levels of actin regulators, such as Rac1 and Arp2/3, activation of the positive feedback loops, precise spatiotemporal regulation of Rho family GTPases, and coupling of protrusion to substrate adhesion 2,62 -then the overall distribution of persistence length can follow truncated power law ( Supplementary Figures 15 and  16). The general conclusion of this modeling effort is, therefore, that synchronization (or desynchronization) of front-back protrusions/retractions may determine the overall motility pattern. The model indicates that cells in which front and back dynamics are synchronized are expected to perform truncated Lévy walks (as in our experiments with metastatic cells) whereas lack of synchronization should translate into diffusive motion (as in non-metastatic cells and metastatic cells treated with inhibitors). Experimentally, we observe some signatures of "desynchronization" of front-rear protrusion-retractions in MDA-MB-231 cells treated with inhibitors (or their combination) that revert Lévy to diffusive motion (see Supplementary Figure 14 for details and Supplementary Movies 16-18). Based on these results, we propose a hypothesis that Lévy-like movement pattern in metastatic cancer cells on 1D microtracks results from the balance/ competition of persistent protrusion-via synchronized frontrear motions dependent on Arp2/3-and an effective myosin IIdependent mechanism for switching the direction of cell motion.

Discussion
Taken together, the above results establish the presence of Lévylike movement patterns in a range of metastatic cancer cells, and also suggest means of changing this motility phenotype. Metastatic cells not only move "faster and more persistently" 10 than their non-metastatic counterparts, but while doing so, also display a fundamentally distinct path structure. This "predatory" navigation strategy is independent of external chemotactic signals as gradients of various chemokines, EGF, or VEGF are all absent in the microtrack experiments. This work contributes to the growing body of knowledge indicating that Lévy-like movement patterns are present and fundamentally important not only for humans 26,27 and multicellular animals [17][18][19][20][21][22][23][24][25] , but also on cellular 6,28,29 and subcellular levels 16 . Lévy walk movement patterns are generally different from and should not be confused with classical models of cell motions such as persistent and simple correlated random walks (PRW/CRW reviewed by 1 ). The latter are generally characterized by ballistic (i.e., persistent in direction) motions at short time and diffusive at long time scales. In contrast, Lévy walks are superdiffusive over long times thus allowing the searcher to move farther away from the starting point in a shorter time than CRW/PRW strategy would. This being said, it should be noted that recent theoretical analyses and experimental evidence from multicellular animal movements indicate that related movement patterns, specifically multi-phasic walks (e.g., composite CRWs and composite Brownian walks), in which the mover switches between two or more kinds of simple walk patterns (mathematically characterized  70 , of the truncated power law exponents μ for cancer cell movements in vivo for the data described in Fig. 3 and 4 (± values give the 95% confidence intervals). Lower cutoff value was a = 5 min for all in vivo data sets. Away from the tumor, in zones 2-4, Akaike weights heavily favor truncated power law with 2 < μ < 3 by step-sizes from two or more exponential distributions which can sum up to a power law distribution 60 ) have been suggested "to have parameters that are fine-tuned to [optimal] Lévy walk" 17,59 . While bi-phasic (or bi-modal) walks have realistic physical correlates (such as run phases and reorientation/wait phases) 64 , it seems difficult to find physical justification for fitting/using more complex multi-exponential distributions to describe step-size distributions, such as the ones we observe here for metastatic cells. In addition, discrete nature of cell motility data we collected precluded reliable fitting of biexponential distributions (see Supplementary Note 2 for more information). Line width = 20 μm, same for all three images. In displacement plots, ten representative trajectories for each treatment are shown. Quantification of motility characteristics (exponents μ, α along with the ± 95% confidence intervals) for MDA-MD-231 cells moving on microtracks and having individual actin-binding proteins inhibited is summarized in Table 3. e Log-log plots of the cells' mean square displacement versus time, x 2 ¼ t α . The slopes correspond to exponents α; note that inhibition of Arp2/3 with CK666 results in diffusive motion, α~1, while inhibition of Myosin II results in ballistic motion, α~2. f The corresponding cumulative frequency distributions, CFDs, of persistence times. Markers are experimental statistics for persistence times. Solid lines are truncated power law fits with respective μ values shown in Table 3. For the chemical inhibitors data shown corresponds to: 40 μM CK666 (for Arp2/3, yellow crosses), 100 μM NSC23755 (Rac1, green triangles), 10 μM Blebbistatin (Myosin II, red rectangles), and control (black circles). The additional results for all drug and siRNA concentrations tested are shown in Supplementary Figures 10-13. See Supplementary Movies 9-12 and [16][17][18] While our observations of Lévy walk motility pattern for metastatic cells moving through very simple experimental environments in the absence of chemotactic gradients suggest that Lévy walk is an intrinsic property of metastatic cancer cells, current work does not exclude the possibility that constrained complex geometries, such as aligned collagen fibers and linear micro-tunnels/tracks, contribute to the emergence of the Lévy walk pattern in complex in vivo environments. In this context, it is interesting that recent theoretical studies have shown that linear constraints, such as one-dimensional micro-channels, provide one of the simplest systems for the emergence of optimal Lévy walks 23,32 . For example, in the system described by Reynolds et al. 23 , a Lévy-like displacement patterns (socalled Weierstrassian Lévy walk 23,29 ) can emerge as a result of the moving object bouncing chaotically off the walls of the channel. While intriguing, such minimalistic mechanical explanation might not be at work in our 1D microtrack system simply because there are no "walls", and cellular motions are determined mostly by cell-substrate adhesive interactions and are gently guided by adhesive micropatterns. At the same time, cellular interactions with complex geometric/mechanical constraints could play role in wider micro-channels/tracks, such as in PDMS micro-channels used by others 65 and in vivo settings 36,43 .
To the best of our knowledge, the work presented here is the first demonstration of Lévy-like movement patterns in adherent mammalian cells, as well as in the context of metastatic cancer. Regarding the latter aspect, we ponder whether adopting the Lévy-like modality of migration might endow metastatic cells with a successful strategy for dispersing and searching for suitable loci where to seed new metastases. The prevailing theory (though disputed by others ,34 ) stipulates that Lévy walk searches are considered to be optimal within a narrow range of circumstances [30][31][32][33] . Specifically, non-destructive Lévy walk searches with μ~2 were considered optimal provided that targets are randomly distributed and scarce and searcher does not have prior knowledge about the target locations [30][31][32] . However, more recent theoretical work has shown that Lévy walk searches (again, with μ~2) are optimal under much broader environmental conditions than previously thought and that Lévy searchers experience fewer long periods of starvation (compared with exponential, composite Brownian and ballistic searchers) 33,66 . Furthermore, the so-called adaptive Lévy searching (where searcher has some knowledge of the target distribution and, upon target detection, responds by switching from extensive Lévy searching to intensive Brownian searching) has been shown to outperform adaptive ballistic and composite Brownian searches 67 . Interestingly, the theoretical 2D path structure analysis 33,34 highlights a remarkable difference in movement patterns between Lévy (with μ~2.5, close to values observed in some of our experiments) and exponential searchers. The former takes many small steps in a focused area resulting in thorough focused exploration (high oversampling) until one large step takes it to another area (similar to 1D path structure of our metastatic cells, see Supplementary . In contrast, an exponential searcher diffuses more gradually and covers the area more evenly without a thorough exploration of any particular spots 33 . This otherwise suboptimal Lévy searcher (μ~2.5) performs extremely well at intermediate levels of prey/target abundance. In the context of metastasis, focused local exploration characteristic of Lévy walk (with μ~2.5) may be useful for finding "hot spots" of defective basement membranes 68 followed by ballistic motions along linear collagen fibers and/or micro-channels toward another set of basement membranes before entering the bloodstream. Interestingly, risk of predation-condition relevant to metastatic cancer cells which during dissemination in vivo are constantly targeted by immune cells-alters optimal search strategy such that searchers executing Lévy walk with 2 < μ < 3 have higher fitness than those with μ ≤ 2 69 . Taken together, our work suggests that Lévy-like movements represent an emergent strategy of how highly motile metastatic cells move within complex microenvironments, and because such movement patterns could be advantageous, they may be maintained and/or fine-tuned by selection pressures.
Naturally, the present results should be generalized to as many types of cells/cancers and conditions (for example, micro-channels) as possible, though especially the in vivo studies are currently limited by the inability to visualize tumors lying deeper into the animal and to obtain long-enough trajectories. In addition, we do not yet fully understand the origin of Lévy walks at the level of cytoskeletal dynamics. While we have developed a simple model (see Supplementary Note 5 and Supplementary Figures 15, 16) that could ascribe these walks by the synchronization of the cell's "front" protrusion (driven by Arp2/3-mediated nucleation of actin filaments) and "back" retraction (enabled by actomyosin contraction), such assumptions need further scrutiny and experimental validation. We see further effort in this area justified by the hope that by eliminating the ability of metastatic cells to execute "predatory" Lévy walks through the human body, it might be possible to decrease the ability of these cells to seed metastatic cancers.

Methods
In vivo observation of cell trajectories. The movements of stable Histone2B/EGF or Histone2B/mCherry expressing non-metastatic B16-F0 and metastatic B16-F10 mouse melanoma cells were observed in live tumors established in mouse dermis by using the dorsal skin-fold chamber model 43 . Dorsal skin-fold chambers were transplanted on a skin flap of C57/B16 J mice (Charles River) and B16-F0 or B16-F10 Histone2B-EGFP/mCherry cells were implanted by injection of 5 × 10 4 -2 × 10 5 cells into the dermis adjacent to the deep dermal vascular plexus. For visualization of blood vessels, AlexaFluor 750-labeled 70 kD dextran was injected intravenously (2 mg/mouse). A customized multiphoton microscope (TriMScope-II, LaVisionBioTec) setup was used which allowed for simultaneous secondharmonic generation (SHG) to reconstruct tissue interfaces, and fluorescence imaging to track blood vessels and the movements of B16-F0 and B16-F10 cells stably expressing Histone-2B/GFP in live mouse dermis. Time-lapse acquisition was done 3-9 days after tumor cell implantation for 1-3 h at 5 or 10 min intervals.
All animal experiments were approved by the ethical committee on animal experiments and performed in the Central Animal Laboratory of the Radboud University, Nijmegen, in accordance with the Dutch Animal Experimentation Act and the European FELASA protocol (www.felasa.eu/guidelines.php).
Data analysis. Heavy-tailed models and fitting procedures: The equations used in the comparisons of the heavy-tailed models can be found in Supplementary  Table 1. Power law and exponential distribution parameters were estimated using analytical expression for maximum likelihood estimator, as discussed elsewhere 56,57 . Parameters for truncated power law, log-normal, and stretched exponential distributions were found using numerical maximum likelihood estimation and verified by visual inspection of likelihood maps. The appearance of largely negative μ values in case of log-normal model is demonstrated in Supplementary Figure 6. Effectively, the log-normal curve is being stretched to better fit power law-like distribution of the data, suggesting that log-normal is not the optimal model for the data. When using only positive μ, the resulting log-normal likelihood and, therefore, wAIC are worse. The resulting parameters are shown in the Supplementary Table 2. Parameter a is a lower cutoff parameter, the choice of which is described below, and b is an upper cutoff parameter used in truncated power law estimation, and corresponds to the maximal point in each dataset.
Lower cutoff estimation: To estimate the choice of the lower cutoff parameter a, we used a technique based on reweighted Kolmogorov-Smirnov (rKS) statistic, described in detail in ref. 56 . In brief, the rKS is calculated for all a values taken from the set of unique values of persistence times in each dataset. Then, a is chosen where rKS is minimal and the number of remaining data points after thresholding by a is not less than 50% of the original dataset. In cases of truncated and regular power law, exponent values in Table 1 are shown for lower cutoffs determined separately using rKS, whereas wAIC for all distributions are calculated using cutoff values for the truncated power law.
For all other experimental details, please see Supplementary Methods/ Supplementary Notes 1-3 online.