Comparative mapping of crawling-cell morphodynamics in deep learning-based feature space

Navigation of fast migrating cells such as amoeba Dictyostelium and immune cells are tightly associated with their morphologies that range from steady polarized forms that support high directionality to those more complex and variable when making frequent turns. Model simulations are essential for quantitative understanding of these features and their origins, however systematic comparisons with real data are underdeveloped. Here, by employing deep-learning-based feature extraction combined with phase-field modeling framework, we show that a low dimensional feature space for 2D migrating cell morphologies obtained from the shape stereotype of keratocytes, Dictyostelium and neutrophils can be fully mapped by an interlinked signaling network of cell-polarization and protrusion dynamics. Our analysis links the data-driven shape analysis to the underlying causalities by identifying key parameters critical for migratory morphologies both normal and aberrant under genetic and pharmacological perturbations. The results underscore the importance of deciphering self-organizing states and their interplay when characterizing morphological phenotypes.


Introduction
Cell migration is a fundamental cellular process that underlies embryonic development, wound healing, immunological surveillance and cancer metastasis. In particular, fast migrating cells such as Dictyostelium and migrating immune cells are versatile in their patterns of movement that range from random exploratory movements with frequent turns to more persistent migration in a straight path. Dictyostelium cells exhibit both random migration [1] as a phagocyte, and more persistent migration when forming a fruiting body. Exploratory interstitial migration in leukocytes [2,3] underlies antigen search and immune surveillance [4], and some are also known to move in a straight line [5]. Frequency of cell turning and their angles is dictated by when and where branched networks of F-actin that drives formation of lateral protrusions called pseudopods occur. In Dictyostelium, pseudopods appear arm-like, and their formation and splitting randomizes cell orientation [6]. Selective maintenance of pseudopods thus provides directional bias in the shallow attractant gradients [7]. Similar F-actin enriched projections in immune cells vary in their appearance from those that are finger-like in DC cells to those more lamellar in neutrophils, however their role in directional choice appears to be conserved [8,9]. On the other hand, ability to move in a straight and persistent manner requires cell polarity which refers to a long-term state having a dominant leading edge enriched in branched F-actin meshwork and a trailing end with crosslinked actomyosin. In certain cells under geometrical confinement, buildup of hydrostatic pressure by contractility can rapidly switch the protrusion to a bleb which is devoid of F-actin [10]. Besides such cases, F-actin driven leading edge protrusion and rear contraction are concomitant in the polarized cells. While the particular shape that cells take depends on the extracellular conditions such as cell-substrate adhesion and diffusible attractants, shapes with broken-symmetry emerge in the absence of extracellular asymmetries, and thus their origins are cell-intrinsic by nature [6,7,11]. The fact that movement of fast migrating cells depend highly on self-deformation contrasts highly to those of mesenchymal cells such as fibroblasts, which move at an order of magnitude slower speed and are strongly dictated by the asymmetries introduced by the adhesive foci.
The common and recurring shapes observed under highly divergent culture conditions and across evolutionary distant species and taxa [12] suggest generality of the self-deforming dynamics in fast migrating cells. Transmigrating neutrophils, genetically or pharmacologically perturbed Dictyostelium [13,14] and certain cancer cells [15] take a canoe-like polarized morphology similar to fish-keratocytes and some protozoan amoebae. Conversely, polarized neutrophils under certain genetic and pharmacological perturbations are known to exhibit increased number of pseudopods [16,17]. These common and interconvertible morphologies suggest that they reflect basic self-organizing states of motile cells that can be quantified and compared with minimal reference to the details of the molecular underpinnings and the permissive extracellular conditions. Due to compounding levels of complexity, quantitative characterization of these canonical morphologies also requires one to leave aside fine-scale protrusions such as filopodia and endocytic cups and apply an appropriate coarse-grained description at the cellular-level. The aspect ratio of fish keratocytes was identified as the major variation in shape features by principle component analysis (PCA) [18]. Variation in more complex features requires other non-trivial measures of characterization. Fourier and related spectral analysis allows one to extract the periodicity in the protrusion-retraction cycle as well as in their spatial ordering [19]. Combined with PCA, Fourier description of Dictyostelium cell shape has shown that the morphologies observed at various steepness of a chemo-attractant gradient can be characterized in a two-dimensional feature space that represents differences in degree of elongation, splitting and polarization [20]. Zernike polynomials in combination with PCA has been used to classify invasive cancer morphologies in two-dimensional feature space [21]. Besides Fourier-based analysis, methods such as tracking of local curvature [22] and pseudopods at the cell edge [7,23] have been employed to characterize spatio-temporal dynamics of membrane protrusions. In RNAi screen of Drosophila culture cells, a large body of handpicked morphology features has been employed to train a classifier by shallow neural networks [24] and Support Vector Machine [25]. These studies indicate that the states of physically realizable morphologies are confined to a relatively low-dimensional feature space [25]. The downside of data-driven approaches, however, are that the analysis often remains in a blackbox making it difficult for one to understand data with reference to the underlying causalities.
A great challenge remains as to how one can quantitatively relate the characterized shapes to the underlying dynamics and vice versa [12,26,27]. For the most basic analysis, it is instructive to formulate a top-down model for isolated cells that is free of extracellular context [12], as behaviors under complex environments may later be deduced, given the repertoire of realizable dynamics, from spatial asymmetries and constraints in the key parameters. In the "graded radial extension model", a polarized morphology similar to that observed in fish keratocytes and neutrophils is described without reference to the underlying mechanism by assuming that the plasma membrane extends radially and that its magnitude is spatially graded along the anterior-posterior axis [28]. Such a steady and graded distribution is thought to result from reaction-diffusion based symmetry breaking in the activity of the polarity signals GTPases Rac and RhoA at the plasma membrane that specifies the state of F-actin at each given place and time. Resource limitation that prevents one state from dominating the other is expected when the sum of the inactive and active form of the small GTPases is approximately fixed in time [29]. Bi-stable reaction-diffusion systems with the above constraint are known to support a protrusive membrane region (front) and contractile membrane region (rear) to co-exist in a spatially separate domains within a cell-a mathematical manifestation of a stable polarized cell shape [29]. On the other hand, pseudopods are transient structures regulated by locally amplified formation of branched F-actin networks. In Dictyostelium this is governed by transient activation in Ras/Rap and PI3K [30], and in case of neutrophils, by Cdc42 and PI3K [31,32]. Because the localized protrusive dynamics occur under uniform conditions, they are thought to arise by noise amplification by excitable regulatory network [23,30,32,33]. Cdc42 in neutrophils and Ras/Rap in Dictyostelium are also known to act positively to strengthen Rac and Rho and hence cell polarity [16,[34][35][36].
Recent mathematical models addressed how interconvertible morphologies of random and persistent migrating cells can be described in a single framework. Modified models of excitability [14,[37][38][39][40] introduce means to prolong propagation of wave-like activities which gives rise to elongated forms with directional persistence. Conversely, in bistable-based models [41], a recent modification to include large noise with memory [42] has demonstrated pseudopodlike dynamics. Despite added complexities, these models still fall within the realm of what excitability or bistability can describe. To which extent these two disparate schemes capture real cell morphologies have not been systematically and quantitatively addressed. In this work, we develop a framework that employs deep learning based classifier to obtain objective measures for shape comparison. By introducing a conceptual model that describes coupling between excitable and bistable regulation, we show that their interplay successfully maps experimentally observed morphologies across the full range of the feature space including those where the existing models fail. Furthermore, the model highlights key parameters that define morphologies under normal and aberrant conditions. The present approach provides a general and extendable framework to characterize varieties of other cell shapes in a data-driven manner which can then be interpreted and tested to further improve hypothesis-driven modeling. Provided that there are large sets of simulated timeseries and real data for feature extraction, the ability to help infer the migratory dynamics from snapshot images should also have practical single-cell applications for cell identification.

A feature space related to cell polarity and pseudopod dynamics can be obtained from classification of stereotype morphologies by deep convolutional neural networks
For systematic extraction of cell morphology features from microscopy data, a convolutional neural network (Fig 1A, lower panel) was trained to classify snapshot mask images of: Dictyostelium (aggregation-stage; 'agg'), neutrophil-like HL-60 and fish keratocytes (Fig 1A, upper panel; S1 Movie). The choice of the reference data was based on the fact that they are well-studied systems and that each represented a stereotype morphology that can be interpreted to represent different degree of pseudopod formation and cell polarity [6,12]. Under our experimental conditions, Dictyostelium cells showed an elongated form in the anterior to posterior direction with locally appearing pseudopods. Fish keratocyte took a canoe-like shape characterized by its long axis orthogonal to the moving direction. HL-60 exhibited an intermediate form between the two where, compared to Dictyostelium, transient protrusions appear less frequently and the overall shape was more horizontally elongated but to a lesser extent than the keratocyte. Image masks of these isolated single cells (Table A in S1 Text) were normalized in size and orientation ( Fig 1A; see Methods). Hyper-parameters for deep-learning were chosen for relative high-accuracy for various network structures. The extracted features were well trained as judged by the high validation accuracy; 97.9% and 89.7%, for the training and the validation data respectively (S1A and S1B Fig; the mean of the last 10 epochs; see Table B in S1 Text for accuracy per dataset). The classification accuracy of the validation data of Dictyostelium, HL-60, and keratocyte is 94.6%, 96.0%, and 87.8%, respectively (Table B in S1 Text). The three nodes F = (F 1 , F 2 , F 3 ) that constituted the second to last layer of the network showed good representation of the three data classes: Dictyoselium (S1C Fig; high F 1 ), HL-60 (S1D Fig; high F 2 ) and keratocyte (S1E Fig; high F 3 ), which were further reduced to two by principal component analysis (PCA). The latency values of PC1, PC2 and PC3 were approximately 66.3%, 33.5% and 0.3%, respectively. The contribution of F 1 , F 2 and F 3 to PC1 = (-0.42,-0.35, 0.84)•F and PC2 = (0.76,-0.64, 0.11)•F indicate that PC1 highlights keratocytelike shapes while penalizing features common to Dictyostelium and HL-60, and PC2 favors Dictyostelium-like features not found in HL-60. Fig 1B shows good separation of the three datasets in the PCA space (~92.8% accuracy) compared to classification based on hand-crafted features (S2D Fig, Tables Q and R in S1 Text). The keratocyte and Dictyostelium (agg) dataset were found in the high PC1 and high PC2 regions, respectively. The HL-60 dataset were mapped to a low PC1 low PC2 region. The feature metrics acquired above, by the very fact that they constitute a good classifier, should be useful to quantify similarity of cell morphology of one's interest in reference to the trained data. The overall relationship between representative cell contours and the morphology feature is shown in Fig 1C. At first glance, higher PC1 appears to indicate more pronounced elongation in the lateral direction while higher PC2 indicates marked longitudinal elongation and protruding edges. These variations in the feature space not only reflected the average morphology differences between the three training classes but also shape changes in time. Fig 1D-1l show representative single-cell timeseries from the validation dataset; i.e. a reserved dataset not used for training. Dictyostelium data showed large fluctuations in both the PC1 and PC2 direction (Fig 1D, 1E, and 1F), whereas the HL60 (Fig 1G, 1H, and 1I) and keratocyte data (Fig 1J, 1K, and 1L) exhibited more marked changes in the PC1 direction (see also S1C-S1E Fig for changes in F). To clarify the relationship between the PC1-PC2 and the extent of elongation, we re-visited the principle component of the fish keratocyte shape variation reported earlier ('shape mode 1' in [18]). Low aspect-ratio shapes (-2σ in 'shape mode 1' [18]) were located near the HL-60 dataset, and high aspect-ratio shapes (2σ in 'shape mode 1' [18]) were located near the keratocyte dataset. Consequently, the shape along the 'shape mode 1' of fish keratocyte [18] constitutes a well-confined manifold in the PC1-PC2 space (Fig 1C; asterisks) to which our fish keratocyte data were also mapped. Furthermore, ellipsoidal shapes with various aspect ratios indicate that an increase in the lateral elongation maps them on the same manifold as the keratocyte dataset ( Fig 1M, PC1 > -1; magenta to pink ellipsoids) while an increase in the degree of head-to-tail elongation maps the ellipsoids to another manifold along the PC2 axis ( Fig 1M,  While the analysis of the stretched ellipsoids clearly demonstrated how the cell orientation and the degree of elongation are mapped to the feature space, the Dictyostelium (agg) dataset was markedly offset from these ellipsoids towards negative PC1 values. To further clarify the nature of this deviation, well-defined polygons were subjected to the same feature analysis (Materials and Methods). All 33 geometrical objects tested were mapped within the region of the PC1-PC2 space spanned by the microscopy data ( Fig 1N). Regular polygons and circles were mapped to a low PC1-PC2 region and their vertically stretched counterparts were mapped to higher PC2 region (Fig 1N, red). Polygons and ellipsoids that were stretched in the lateral direction marked high PC1 value (Fig 1N, blue; S1H and S1I Figs). Of particular note were the star objects which scored highest in the PC2 value and deviated markedly in the feature space from the ellipsoids and other objects when stretched (S1I Fig). Stars have higher PC2 compared to squares and triangles indicating that PC2 reflects pointed edges. Comparisons between upright and vertically flipped stars and triangles indicated that degree of pointedness towards the cell front also affects PC2 but in a complex way (Fig 1N; see also S1J Fig). An analysis of more asymmetrical geometries with varying number of edges showed that they map to a domain in high PC2 with high variations towards negative PC1 values as in the Dictyostelium (agg) dataset; i.e. away from ellipsoidal shapes (Fig 1O). Rotating the multi-edge network (lower panel). Masks were normalized in area and aligned downwards in the migrating direction. The feature vector F = (F 1 , F 2 , F 3 ) defined by the three nodes in the last layer was further reduced to PC1 and PC2 by PCA. (B) Mapping of trained data in PC1 and PC2 (dark green + (aggregation-stage Dictyostelium), dark red + (HL-60) and yellow + (keratocyte)). Each data point represents timeaveraged scores from a time-series of a single cell. (C) Representative cell contours mapped to the feature space. Asterisks are the first principal shape variation of fish keratocyte (images taken from [18] (-2σ (blue) to +2σ from the mean (purple)). forms with the small PC1 values (Fig 1O; 0 to 30 degrees) bring about decrease in PC1 (S1K Fig). Further rotation increases PC1 (S1K Fig; 60 to 90 degrees) due to the shape now appearing more horizontally elongated overall. Generality of the results was confirmed in independent real-cell data by analyzing fully differentiated prespore cells of Dictyostelium which is elongated longitudinally and lacks pseudopods thus mapping identically to the elliposoids (S3A and S3B Fig). Likewise, effector T cells with their signature branching protrusions were broadly distributed along PC2 (S3C Fig, Th1) compared to markedly less polarized regulatory T cells (S3D Fig, Treg). These analysis indicate that our classifier was able to yield data-driven representation of complex signatures with respect to the cell orientation for both local pseudopodal protrusions and more global cell elongation.

The 'ideal cell' model recapitulates a generalized morphological landscape constrained by the choice of the protrusion speed and the balance between the local protrusion and global polarity
Let us now introduce an ideal cell model (see Equations: Eqs 1 and 2) that serves as a canonical shape generator of migrating cells that is strictly constrained by the transient protrusion dynamics and polarization (S1 Text). To describe the interfacial membrane mechanics, we employed the phase-field method with the addition of an active force F prot = a w W (Eq 1; Table C in S1 Text). The model also consists of spatio-temporal dynamics of variable W and variables U, V that define global cell polarity and local protrusions, respectively (Fig 2A; Eq 2). These variables are abstract representation of how the respective signaling molecules-such as Rac, Rho for W and Ras, Cdc42, PI3K for V-are regulated in space and time (Materials and Methods). W is bistable and takes either a low or a high state which signifies the retracting rear and the expanding front, respectively. The low state indicates that W has converted to the other form W � while the integrated sum of W and W � is fixed to W tot . Owing to this constraint, the dynamics of W exhibits wave-pinning and therefore supports a polarized shape with W being high at one end and low at the other end ( Fig 2B, t = 660s). In addition, there is an excitable network that describes conversion of U to V which is invoked by small signal fluctuations. The important assumption in the model is that these two core networks are coupled. On top of the core bistability in W, V promotes conversion of W from the low to the high state, and W catalyzes amplification of V (Fig 2A). Because V amplifies noisy fluctuations and generates local protrusions through W, a leading edge defined by a region with high W is most likely to be perturbed and often split into two ( Fig 2B; t = 700). This is however transient, as W by itself works to maintain global unipolarity hence only one protrusion survives ( Fig 2B; t = 800s). In some cases, a new protrusion can also form away from the anterior and more towards the lateral side and still develop into a new dominant front ( Fig 2C). These features required full 3-variable equations (Eq 2; Table D in S1 Text) and were of particular importance in our comparative analysis. Insofar as our parameter search (Tables E and F [37,43,44].
In total, 228 parameter conditions (Table G in S1 Text) were selected by heuristic sampling where we performed grid search in the parameter subspace around hand-picked reference points (S5 Fig; Material & Methods). The advantage of the 'ideal cell' model over the existing models focused on specific cells and conditions is that, it can describe a generalized Schematics of the dynamical model. Self-amplifying and excitatory synthesis from U to V [70] induces protrusion factor W which can be prolonged by the bistable dynamics (Eq 2). The reaction takes place in a region ϕ = 1 governed by interface physics according to the phase-field equation (Eq 1). The membrane expands outward into an unoccupied region ϕ = 0 in the direction perpendicular to the border at a rate proportional to local W. Main parameters (χ U , γ, μ, k W1 , ρ, a W ) are denoted along the associated reaction steps. (B, C) Representative model behavior (overlay plots for U, V, W); front splitting (B) and pseudopod formation morphological landscape incurred by the choice of the protrusion speed and the balance between the local protrusion and global polarity. As we shall describe later, this versatility allows us to resolve and characterize cells in different developmental stages and under perturbations. First, of particular interest was whether the cell was elongated in the lateral or anterior-posterior direction. We found W tot has a large influence in the direction and the degree of overall cell elongation ( Fig 2D; see also S2 Movie). For small W tot , the system behaves as a passive interface that relaxes to the state of minimum surface energy of fixed area, thus the simulated morphology was near circular ( , where the simulated morphology was again near circular because W occupied an entire cell. Accordingly, moderately high W tot (W tot = 80) supported relatively high F 1 and high F 2 on average (indicating shapes resembling Dictyostelium and HL-60) ( Fig  2D; bottom panel). In the PC1-PC2 space, small W tot yielded low PC1 and low PC2, whereas at a moderately high W tot , PC2 took high values (Fig 2E). At high W tot , the average PC1 and PC2 decreased. The overall dependency on W tot were conserved when other parameters were varied ( Fig 2E; magenta and blue). The other important parameters that affected cell polarity was ρ; the strength of autocatalysis in the interconversion reaction between W and W � . High ρ means that the non-zero roots of the cubic equation −ρW 3 +ρW � W 2 −W = 0 is large thus supporting a larger domain with high W. Therefore, at low ρ, the leading edge is small and cells become elongated in the moving direction (S5A Fig; ρ  Occurrence of local protrusions depended largely on k W1 and a W . k W1 specifies the depth of the bistable well. Because k W1 by definition defines how strong protrusive activity is dominated by bistability of W rather than through V, it can be viewed as representation of relative activity in Arp2/3 and crosslinker Myosin II. For small k W1 , the front-rear asymmetry was weak, and the overall cell shape was near circular (Fig 2F; k W1 = 1.1). Intermediate k W1 gave rise to mixed dynamics where local protrusions induced asymmetrical deformation however without persistent front-to-back polarity ( Fig 2F; k W1 = 2.2-55). Large k W1 elevates W which makes it less affected by the dynamics of U and V, and thus supports elongated shape with more marked polar asymmetry in W (Fig 2F, k W1 = 110; S3 Movie). Accordingly, we obtained high F 1 (i.e. high resemblance to Dictyostelium (agg)) ( Fig 2F, lower panel), and hence high PC2 (Fig 2G). Deletion of myosin light chain reduces pseudopodia and strengthens polarity in Dictyostelium [45]. Conversely, increased myosin light chain kinase expression is known to reduce lamellipodia size and induce multiple protrusions in keratocyte independently of Rho kinase and membrane tension [46]. The appearance of local protrusions also depended strongly on the protrusion force a W (

PLOS COMPUTATIONAL BIOLOGY
Model-based analysis of crawling-cell morphology in deep learning-based feature space cell orientation changed frequently (S4 Movie). There was an increase and a decrease in F 1 and F 2 respectively (i.e. high resemblance to Dictyostelium (agg)) ( Fig 2H, lower panel). The PC2 score increased accordingly (Fig 2I).
Parameters that affected the pseudopod dynamics were μ and γ which define the downregulation rate of V and U, respectively. Excitable dynamics in neutrophils and Dictyostelium are associated with PIP3, and thus μ can be viewed as activity of phosphatases and kinases other than those directly involved in conversion between U and V. Low μ (Fig 2J; μ = 0.1; S5 Movie) elevates V, hence the concomitant increase in W supported a laterally elongated shape. Dictyostelium cells are known to take similarly elongated form when perturbed with 5-phosphatase that effectively reduces plasma membrane PTEN and increases PIP3 [14]. The polarized shape was highly persistent as high V renders the patterning less prone to noise perturbation. On the other hand, at intermediate to high value of μ, cells became more elongated longitudinally and the polarity was less persistent (Fig 2J; μ = 0.5-0.9). Here, the high W domain was easily disrupted; fronts frequently split, and the cell orientation was altered (e.g., a Y-shaped front in Fig  2J at μ = 0.7 and 0.9). Accordingly, F 3 (Fig 2J; bottom panel) and the PC1 score decreased at high μ (Fig 2K). Similarly, at low γ, pseudopods split frequently and new pseudopods were rare

Morphology-based mapping of model parameters can help infer candidate dynamics
The distribution of the simulation data in the PC1-PC2 space spanned a large region occupied by the training dataset ( Fig 3A, black circles) further vindicated the ability of the ideal cell model to describe the characteristic morphologies. Proximity of the time-averaged simulated morphologies ( Fig 3B) to the average of the three reference dataset was analyzed by computing the Euclidean distance in the feature space F = (F 1 , F 2 , F 3 ). According to the reference data, the distance was designated as Score-D (Dictyostelium), Score-H (HL-60), Score-K (keratocytes), and ranked in the ascending order; i.e. a low score means high similarity (Table H in S1 Text). The time averaged morphology feature in the PC1-PC2 space and the time-series of the top ranking simulations are shown in Fig 3B (filled circles) and Fig 3C-3K, respectively. Simulations with high Score-D on average exhibited morphology that closely resembled the aggregation-stage Dictyostelium with their elongated form in the anterior-posterior direction accompanied by a few pseudopods that frequently reoriented cell directionality (Fig 3C; S6 Movie). Similarly, simulations that ranked high for Score-H (Fig 3F; S7 Movie) exhibited fanlike cell shape that moved directionally with some occasional turning as observed in HL-60. These high ranking parameter sets were found near the median of the reference dataset in PC1-PC2 space (Fig 3E, 3H, and 3K). When high-ranking simulations were retested for classification per snapshot, close to or higher than 90% were correctly classified as Dictyostelium (Table I in S1 Text) or HL-60-like (Table J in S1 Text). As for high Score-K simulations, the simulated shapes were canoe-like with high directional persistence (Fig 3I; S8 Movie). On average, they deviated transiently in the PC1 direction thus making classification per snapshot less clear (< 66% accuracy). The deviation was due to occasional shape perturbation by  Table M in S1 Text).
Although similarity was evaluated based on still images, dynamics of high ranking simulations were by and large consistent with those of real cells. This was well illustrated in the kymographs of boundary curvatures (Fig 4A and 4B upper panels); the anterior and posterior regions can be identified by high curvature regions with either positive or negative velocity respectively (Fig 4A and 4B lower panels). Both Dictyostelium and HL60 (Fig 4A and 4B left panels) showed anterior projections (Fig 4A and 4B bottom panel red regions) that bifurcated from time to time and traveled towards the high curvature region at the rear [47]. The wavelike appearance was somewhat more prominent in Dictyostelium. In both cell types, the posterior end was characterized by a high curvature region that persisted over time. These dynamical features were well recapitulated in the simulations (Fig 4A and 4B right panels). Moreover, there was good agreement between the simulations and the real data in the cell trajectories. The mean square displacement (MSD) of the centroid showed a characteristic time-scale dependency where it was proportional to the square of the elapsed time (ΔT 2 ) for ΔT < τ 0 ( Fig  4C, magenta line) and to ΔT for ΔT > τ 0 (time domain, Fig 4C red line). In other words, cells moved ballistically i.e. at a constant velocity for ΔT < τ 0 and more like a Brownian particle for ΔT > τ 0 . τ 0 can be interpreted as the persistence time for directional migration, and square root of the MSD at the inflection point X 0 characterizes the persistence length. Throughout this paper, we chose the time-scale factor τ 0 = 10 based on approximate matching in the crossover point of the two regression lines between the top ranking simulations and the real-cell data (Fig 4C and 4D; red lines). For the top Score-D simulations, we obtained τ 0 = 87 sec, X 0 = 9.8 cell length, compared to τ 0 = 151 sec and X 0 = 15.2 cell-length in the real data which are in good agreement with values reported earlier [48]. Trajectories of top ranking simulations for Score-H were more persistent (τ 0 = 257 sec, X 0 = 22.1 cell length) as was the case for the real HL60 data (τ 0 = 278 sec, X 0 = 47.7 cell length).
The other important feature of random migration is the relation between pseudopod dynamics and the cell orientation [7]. New pseudopods frequently appeared in vacant regions (Fig 4E left), or on top of a pre-existing pseudopod thereby giving the cell the appearance of Yshape ('Y-split' in [7]). In other cases, pseudopods continued to extend while turning ('oneway-split' in [7]). The relative occurrence of de novo formation of pseudopods was approximately 47% for Dictyostelium (agg) data and 34% for the high ranking simulations (Fig 4F). For HL60, de novo formation in real cell data was 12% and in the high ranking simulations 21% (Fig 4F). The extension angles relative to the direction of centroid displacement was about 20-40˚for both Dictyostelium and HL60 data and their top-ranking simulation counterparts (Fig 4G and 4H). Although there was some overrepresentation of extension angles around 90˚in the simulation, the angles above 120˚were rare in both real data and simulations. All in all, these results demonstrate that the model, albeit its simplification, is able to recapitulate semi-quantitatively both the persistent random walk behavior and the underlying morphology dynamics in Dictyostelium and HL-60 cells.

Mapping of morphological diversification predicts key parameters for state transition
Although the datasets analyzed above showed little overlap with one another in the feature space, it should be noted that these coordinates are by no means singular representation of specific cell-types and species from which the data were obtained. As we saw above, there was a large cell-cell variability in the fish keratocyte data that constituted a distinct manifold in the feature space (Fig 1B; yellow). Likewise, cell-cell variability was evident in the aggregationstage Dictyostelium cells along the PC2 axis (Fig 1B; green). To see how changes in cell-intrinsic properties alter their positions in the feature space, data from new experimental conditions expected to alter cell polarity were studied (Fig 5A). Undifferentiated (vegetative) Dictyostelium cells took less elongated shape than the aggregation-stage Dictyostelium cells under the same substrate and buffer condition. Their aspect ratio on average was smaller than aggregation-stage Dictyostelium but larger than that of HL-60 ( Fig 5B). Accordingly, in the PC1-PC2 space, the vegetative Dictyostelium was mapped between aggregation-stage Dictyostelium and HL-60 ( Fig 5A; magenta reverse triangles). Model simulations that ranked similar to the vegetative Dictyostelium data (Fig 5C; S9 Movie) had small a w in common (Table L in S1 Text). Similarly, we analyzed HL-60 cells treated with microtubule destabilizer nocodazole which is known to strengthen neutrophil cell polarity [49,50]. Nocodazole-treated HL-60 cells showed morphology similar to keratocyte with a somewhat smaller aspect ratio ( Fig 5D) and were mapped between the non-treated HL-60 and the keratocyte datasets (Fig 5A; blue triangle). The nocodazole-treated HL-60 cells exhibited shape fluctuations making them wobble which was also observed in the simulations (S10 Movie). These features were well represented in the respective simulations that ranked high for shape similarity (Fig 5E). In addition to high W tot , high ranking simulations had relatively low a w or k w1 (Table L in S1 Text). This can be interpreted from the fact that low a w prevents a cell from breaking apart at high W tot , while low k w1 makes the polarized front less pronounced and more sensitive to noise perturbation.
Next, a null strain of racE (Dicytostelium RhoA homologue) was chosen for the analysis because of its aberrant cell shape that resembled fish keratocytes [51]. Under our experimental condition, relatively undifferentiated racE-cells exhibited canoe-like shapes that were similar to the fish keratocyte but more dynamic (Fig 5F). Small fragmented pieces were observed to occasionally split from high curvature regions (Fig 5F). These data marked relatively high PC1 and PC2 scores and were mapped between the keratocyte and the aggregation-stage Dictyostelium data (Fig 5A; black square). The top ranking simulations for the racE-data (Table L in S1 Text; high Score-D (racE-)) exhibited remarkably similar morphology dynamics characterized by high lateral deformation and occasional fragmentation (Fig 5G; S11 Movie). As expected, there was a large overlap with high Score-K data, leaving only the top 3 Score-K data ( Table H in S1 Text bottom rows; ||F|| 2 < 3 x 10 4 ) that uniquely mapped to the keratocyte data. Compared to either the high Score-D (veg) or the high Score-D (agg) data, the high Score-D (racE-) data had large W tot (Table L in S1 Text) consistent with the laterally extended cell shape. All in all, the above analysis suggests key parameters a w and W tot that are pivotal for the state transition between the characteristic morphologies. This was further verified by studying unexamined regions in the high-dimensional parameter space near the top ranking simulations for vegetative-stage Dictyostelium. Increasing the value of a w brought the morphology score closer to that of the aggregation-stage cell (Fig 5H). In contrast, increasing W tot increased PC1 and brought the shape state closer to the nocodazole-treated HL-60 and racEcells (Fig 5H).

Discussion
Deep learning-based extraction of shape features is a promising avenue for cell classification and identification [52]. In this study, we presented a hybrid cell morphology analysis that combined deep-learning-based feature extraction and dynamical model simulations. Convolutional neural network was trained to classify stereotype migrating cell morphology represented by Dictyostelium, HL-60 and fish keratocytes. The feature vector of the trained classifier showed that the three representative morphologies can be described in low dimensional feature space; the first component represents elongation in the lateral direction and the second component represents anterior-posterior elongation and edges. Earlier studies have addressed a low dimensional embedding of morphology with milder variations [18,19,[20][21][22]24,25]. The present results extend this approach to encompass more divergent morphologies commonly observed in different migratory cell-lines. Our new feature space essentially expanded the principle shape variations associated with cell polarity previously identified in the keratocyte data [18] to encompass narrowly polarized cell morphologies with pseudopodal extensions. The finding of low dimensional feature space is in line with an earlier feature representation of Dictyostelium cells based on Fourier-mode decomposition of the cell contour [20]. While being less analytically clear than the Fourier analysis, the present approach allows one to distinguish morphology with respect to the cell orientation and thus suited to analyze migratory cells. Since it did not take time sequence into consideration, a potential shortcoming of the analysis is that the comparision with real cells had to be amended with other dynamic features such as centroid movement and pseudopod splitting. Another downside is that application to cells with no clear orientation is not straightforward. Interestingly, sample conditions expected to alter cell polarity-vegetative stage Dictyostelium, racE-/Dictyostelium, and the nocodazole treated HL-60 cells were mapped to intermediate coordinates spanned by these two major manifolds. Moreover, the new sample conditions were each found clustered without significant overlaps with the training data sets. These results suggest that the region spanned by the three training datasets (keratocyte, H60, Dictyostelium (agg)) in the PC1-PC2 coordinate constitutes a space continuously occupied by forms realizable by genetic and phenotypic variations.
Search for the core regulatory logic behind cell deformation is important as it provides a much-needed basis to bridge the extracted shape features and the dynamics. Recent models have hinted at how the common shapes of migrating cells and their interconversions may come about (Table S in S1 Text). Excitable systems [14,[37][38][39][40] describe pseudopod-like dynamics however at the expense of poor description of persistence; they exhibit unusually flat and periodically varying morphologies due to the wavefront dynamics (S1 Text; S4A-S4D, S9 and S10 Figs). Bistable models depicts polarity well, however simple addition of excitatory protrusions [42,53] do not recapitulate morphologies associated with formation of new pseudopods (S1 Text; S11 Fig, Tables T and U in S1 Text). These analyses show that a simple additive approach to complement the missing part is insufficient. Our model therefore takes into  Table N

PLOS COMPUTATIONAL BIOLOGY
Model-based analysis of crawling-cell morphology in deep learning-based feature space account the competition and feedback between the excitable protrusion and bistable polarization explicitly. Amplification of noisy signals through a feedback loop consisted of the Wdependent amplification of V (Eq 2C) and amplification of W by V (Eq 2B) supports de novo formation and splitting of pseudopods while maintaining the global cell polarity. These dynamics compete for dominance under limited W tot i.e. the maximal protrusive activity. The resulting morphologies capture features that were not fully represented by earlier models (S1 Text).
Given its abstract nature, the ideal cell model is tailored to describe generic migratory morphologies observed across cell-types and species at the cost of being limited in predicting how exactly the coupling is implemented biochemically which may differ depending on the system. In neutrophils, excitable leading edge dynamics are governed by Cdc42 [32], and it also acts globally to enforce cell polarity by promoting actomyosin contractility through its effector WASP in a microtubule dependent manner [54]. In Dictyostelium, excitable dynamics at the leading edge is observed at the level of Ras which also interacts with GDP-bound form of RacE to strengthen cell polarity [36]. In this respect, mapping of the racE-morphology to a high W tot state (Fig 5H) hints at the nature of the competition at least for Dictyostelium with regard to the states of F-actin: the contractile cortical meshwork that are crosslinked with myosin II or the protrusive dendritic meshwork that requires the Arp2/3 complex for side-branching nucleation. RacE is essential for plasma membrane localization of Diaphanous-related formins (DRFs) [51], and deletion of DRFs (ForA-/ForE-/ForH-) results in the loss of cortical actin. Since the morphology of the null mutant of DRFs phenocopies that of the racE-null cells [51], the increase in W tot are associated with the absence of DRFs from the plasma membrane. On the other hand, the fluctuating protrusions are largely associated with fast idling pulses of Scar/ Wave activities which are amplified by the excitable network [55]. Recent studies suggested that actin nucleators such as formins and Arp2/3 are competing for a limited pool of actin monomers and/or their upstream activators such as Rac-GTP [56][57][58]. Such notion is also supported by an observation that the amount of F-actin is compensated in Scar-/WASP-cells by increased localization of ForH at the cortex [59]. Taken together with our mapping of racEdata, these observations are in line with our current model view that excitability and bistable regulatory networks compete for dominance over limited W tot .
The variations in the distinct morphologies of differentiating Dictyostelium cells suggest alterations in the key parameters that serves as a control point. The difference between the vegetative-and aggregation-stage Dictyostelium was ascribed mainly to an increase in the membrane protrusion force a w . The increase in a w can be understood from the fact that Rac1 [60] and SCAR [61], the essential factors for Arp2/3 activation, are known to be expressed at low levels in the vegetative-stage then increase markedly in the aggregation-stage cells. A recent study based on an excitable model [37] suggested a progressive state transition from a circular to amoeboid then to a keratocyte-like shape by the increase in the protrusion force (Fig 2F. in [37]). Rather, our model predicts that changes in the protrusive force should allow a direct transition from a circular shape (low PC1 low PC2) state to either amoeboid or keratocyte-like form depending on W tot . Such direct transition has been demonstrated experimentally and was attributed to an increased activity of a nested excitable network [14]. However, as discussed above, the elongated shapes simulated in their model were oscillatory and lacked persistency (S10 Fig). The presence of the polarity dynamics is essential for the persistent and longitudinally extended cell shapes. Relatively low D w is required in both vegetative and aggregation-stage Dictyostelium to prevent the polarity dynamics from completely winning over the excitable dynamics. A highly polarized form at high D w ; i.e. the 1-variable model limit z/ Dw ! 0 (S4E- S4H Fig; S1 Text) is indeed reached by cells that further differentiated into prespore cell-type (S4I Fig). While it is possible that certain cells are in a decoupled state (z = 0), the requirement of large D w for cell polarity signifies the importance of W acting globally. Large D w may not necessarily be mediated by pure diffusion as the present model postulates, but instead could be realized by other transport processes implicated in cell polarity such as membrane flow or the myosin-II dependent global actin flow. Global actin flow has been shown to maintain asymmetric distribution of de-filamenting factors [62,63], however such global flow may not always be present in polarized cells [64]. Since these transport processes are tied to cortical actin, they are naturally accompanied by changes in membrane tension [65] which should also be part of the feedback process from W to V.
Quantitative and systematic analysis of model outcomes will only increase its importance as we proceed further to unwind causality behind detailed geometries and dynamics associated with specific cells and conditions. The present framework of data analysis potentially provides means to test and improve specific models of migrating cells. For example, our ideal cell model gave rise to pseudopods from the tail region more frequently compared to the real Dictyostelium (agg) data. The discrepancy could be due to the fact that retraction was assumed to be driven only by the area conservation and that no regulated contractility was explicitly described. While this approximation can be justified when there is reciprocity between the front expansion and the rear contraction as has been shown to hold independently of actomyosin in neutrophils [50], the present model could be modified in the future to include local cortical actomyosin regulation when analyzing detailed shapes of the cell rear and the blebbased front protrusion. Further improvement of the model and increasing dimensionality of the feature space may work hand in hand with extending the present analysis to classify morphologies exhibited by other cell types of wildtype and mutant backgrounds. For example, the present analysis fails to distinguish the pancake-like shape known for Rac and Rap related mutants that result from uniform expansion [44] and similarly round (i.e. low PC1, low PC2) cells inhibited of actin polymerization. This limit maybe overcome by introducing absolute size instead of normalizing the area so as to distinguish spherical cell versus flattened cell in two-dimensional cell masks. Expanding the analysis to 3-dimensional images would also be better suited to the present machine learning approach. As resolutions and dimensions are increased, the cell-shape based analysis may be supplemented with fluorescence image data of cytoskeletons and their regulators. Given the significant bottleneck in the present simulation by the huge computational loads which required parallel computation by GP-GPU, other avenues of coarse-graining maybe required to extend the present approach to a larger multimodal analysis.

Dictyostelium and HL-60 cell culture and data acquisition
Time-lapse data of freely migrating Dictyostelium, neutrophil-like HL-60, fish keratocyte and differentiated T mouse cells were acquired with an inverted microscope using either 20, 40 or 60x objective lens. For Dictyostelium and HL-60, cells expressing Lifeact fused to mNeonGreen [66] and mTurquoise2 were employed, respectively. A Lifeact-mTurquoise2 expression vector was constructed by ligating Lifeact-mTurquoise2 into an episomally replicating plasmid pEB-Multi Neo (WAKO, 057-08131) at restriction sites XhoI and NotI. The Lifeact-mTurquoise2 expressing stable HL60 cell line was obtained by introducing the plasmid by electroporation (NEPA21; Nepa Gene, Ltd., Chiba, Japan) followed by G418 selection (1 mg mL -1 ) after 2 days. For fish keratocytes, DIC images were employed. Dictyostelium cells were grown axenically and obtained according to standard protocols as previously described [67]. Vegetative Dictyostelium AX4 cells (N = 1694 snapshots from 18 timeseries;), aggregation-stage Dictyostelium AX4 cells (starved for 3.5 hours; N = 2841 snapshots 19 timeseries;), vegetative LifeactGFP/racE-cells (N = 330 snapshots from 8 timeseries;) were plated on a non-coated coverglass and images were acquired at 5 sec/frame (aggregation-stage AX4 and racE-cell) or 15 sec/frame (vegetative AX4). Neutrophil-like HL-60 cells were grown in RPMI1640/glutamate media (Wako 189-02145) supplemented with 12% FBS (Sigma 172012). Differentiated HL-60 cells were obtained by treating the cells with DMSO for 3 days. Images of HL-60 cells on fibronectin-coated glass plates in the presence of 1 nM fMLP (N = 3468 snapshots from 23 timeseries;) were taken at 5 sec/frame. For nocodazole treatment, differentiated HL-60 cells were collected by centrifugation, suspended in fresh HBSS containing 20 μM nocodazole and plated on a coverslip pre-coated with 1-2% BSA in PBS. Data were acquired within 20-75 min in the same medium in the presence of nocodazole (N = 181 snapshots from 3 timeseries).

Deep-learning-based feature extraction
Mask images were pre-processed as follows (Fig 1A, upper panel): (i) the migration direction was determined from the centroid displacement at a five timeframe interval (equivalent to 25 sec for aggregation-stage and racE-Dictyostelium data, 10 sec for keratocyte data) except for vegetative Dictyostelium data where 1 timeframe (equivalent to 15 sec) was used, (ii) binarized mask image was rotated to align the migration direction to the y-axis, (iii) the image was rescaled so that the cell area is equal to the area of a circle with 25 pixel diameter. The rescaled masks were each embedded at the center of a blank square frame of 64x64 pixels. The exact spatial resolution of mask images varied from sample to sample due to rescaling, however they were all in the order of~0.5 μm/pixel. Convolutional neural network (Fig 1A, bottom panel) was implemented using Keras (https://keras.io) with TensorFlow backend. To make the sample size of the three datasets near equal, data augmentation was performed by rotating the original masks at angles (< ±5 deg) randomly picked from a uniform distribution (see Table A in S1 Text for the number of samples). Input vectors were processed through layers of convolution operation (Fig 1A, bottom panel; 'convolution') in addition to layers of max pooling operation with a 3x3 kernel to render the analysis robust to positional deviation (Fig 1A, bottom  panel; 'pooling'). These were then processed through a set of densely connected layers with rectified linear and hyperbolic tangent activation function (Fig 1A, bottom panel; 'rel' and 'tanh'). In the final layers, the dimension of the vector was reduced to three and were passed to 'softmax' activation function. The values of the three nodes (F 1 , F 2 , F 3 ) before the final softmax layer were employed to represent cell shapes. The number of training epoch was 2000 which was sufficient for adequate learning as determined by the accuracy and loss values (S1A and S1B Fig). The vector representing three nodes F = (F 1 , F 2 , F 3 ) was further reduced in dimension by principal component analysis (PCA). The PCA parameters were acquired from the average F of 54 timeseries in total (19,23, and 12 samples for Dictyostelium, HL-60 and keratocyte, respectively).

Geometrical analysis of feature space
To examine mapping of 2-dimensional geometries in the feature space, 7 well-defined objects were used; circle, isosceles triangle, asymmetric right triangle, rectangle, rhombus, pentagon and star polygon in a specified orientation (Fig 1N). Triangles, pentagon and star were flipped vertically to obtain total of 11 basic objects. The library was further expanded to 33 shapes by including variants of the basic geometries with different aspect ratio 1:1, 1:2 (Fig 1N red, horizontally elongated shape to the moving direction) and 2:1 (Fig 1N blue, flattened shape with vertically elongated to the moving direction).

Model equations
To numerically simulate the interface between the plasma membrane and the extracellular space, we employed a phase-field equation in the following form [42,69,70].
The equation describes the dynamics of a continuous state variable ϕ(r;t) in a twodimensional space that specifies whether a position r is occupied (ϕ = 1) or not occupied (ϕ = 0) by a cell at time t. In the present study, we consider initial conditions with a single continuous domain with ϕ = 1. The parameter τ is the viscous friction coefficient. The first term of the r.h.s represents effective surface tension, where Δϕ and G are derived from membrane energy. G is Landau functional describing a bi-stable potential. Here we chose G(ϕ) = 18ϕ 2 (1−ϕ) 2 [69]. The second term describes restoring force that keeps the cell area close to A 0 . The Eq 1 assumes that the bending energy is negligible [70]. The third term represents active force with magnitude F prot that is perpendicular to the boundary |rϕ|6 ¼0 and thus drives membrane extension. In the present simulations, parameters in Eq 1 were set so that they are an order of magnitude within generally accepted values (Table C in S1 Text); surface tension η = 1.0 [pN], [71], cell area A 0 = 78.83 [μm 2 ] (~5 μm radius circle), protrusive force by actin polymerization a w W = 0.8-4.0 [pN/μm] (for W~1) [72,73]. Since τ was not well constrained experimentally and expected to differ between cell-types and the culture conditions, we adopted an empirical value 0.83 [pN/μm 2 ] [70], which was then calibrated retrospectively by a multiplier τ' so that the time scale of the simulated cell trajectories match with that of real data as described in the later section. The area constraint M = 0.5 pN/μm 3 and the size of the boundary layer ε = 1.0 [μm] were set close to those in the earlier studies [69,70].
The 'ideal cell' that can take close to all of the basic phases of morphology features that we have examined (Fig 1) should consist of two main features: 1) transient appearance of localized protrusions and 2) prolonged presence of single expanding edge and retracting tail to appear under homogeneous extracellular conditions [1]. Here, we formulate a reaction-diffusion model that describes these two processes mathematically as follows: where the first two equations for U and V describe an excitable reaction network for transient protrusive dynamics, and the third equation for W describes polarization dynamics, respectively (Fig 2A; S1 Text). The role of the excitable reaction network is to generate transient signals for local protrusions in a minute to a few minutes time-scale by amplifying small edge fluctuations of seconds order [55]. In neutrophils, excitable dynamics of Cdc42 and PI3K activity [32,33] is essential for front protrusions. In Dictostelium, excitable dynamics are observed at the level of spontaneous Ras and PI3K activation [30,70,74]. Here, we adopted equations originally introduced to study excitable PI3K activities and the resulting F-actin waves in Dictyostelium (Eqs 2A and 2B) [70]. Parameters α and β are the rate constants of reaction U ! V and V ! U multiplied by the time-scale factor τ 0 . The source of edge fluctuation is introduced as a noise term ϕ N(r,t) [70] in Eq 2B which is amplified through V by a positivefeedback described in the first term. Increase in V is then slowed down due to depletion in U, and the system eventually recovers the original resting state. The expanding membrane region is determined by W governed by Eq 2C (Fig 2A; see S1 Text for derivation) which is similar in form to the wave-pinning [29]. The same type of equation has been used earlier to study polarized cell shape in fish keratocyte [75] and Dictyostelium [42]. The first term describes reaction kinetics with bistability, and the second term describes diffusion. We assume that the sum of W and its reciprocal state W � is conserved (W tot = const.) [29] and that W � diffuses much faster than W and thus can be approximated as uniform in space. Hence we obtain Due to the global constraint, this class of bistable system gives rise to coexistence of low and high W regions due to stalling of a transitory wave front-a property known as "wave pinning" [29,75]. For the sake of brevity, we shall embed rear retraction passively in the form of restoring force. Thus we set F prot = a W W in Eq 1 [69] so that the edge expands where W is high, and the rest of the domain with high W tot /A 0 -W as a result contracts due to area conservation imposed by the second term in Eq 1. Note that W only specifies where the protrusions and retraction take place and does not assume the origin of their driving force. The variable W can thus be interpreted as an abstract representation of bistable signals such associated with the leading edge of polarize cells as Rac-GTP, or alternatively, W tot /A 0 -W with the reciprocal spatial profile to indirectly represent signals that regulate rear contraction such as RhoA-GTP. A recent study has shown that Rac-GTP accumulates at the front of migrating neutrophils regardless of whether it is F-actin filled or blebbing [65]. Based on the inherent ability to maintain persistent front-and-back unipolarity, the polarity network (Eq 2C) serves as a spatio-temporal filter to select or remove local protrusions. This coupling is described by assuming that W positively regulates the positive-feedback amplification of V (Eq 2C; r.h.s first term) in addition to W feeding back to increase V (Eq 2A) to reflect F-actin-or tension-dependent positive regulation of leading edge signals in neutrophils and Dictyostelium [65,76,77].

Parameter search
Model parameters were selected for the systematic feature analysis (Fig 3A; Table G in S1 Text) based on the following criteria: Of the total 22 parameters (Table C and D in S1 Text), eight parameters (W tot , a W , k W1 , μ, γ, ρ, D W , χ U ) were chosen for detail analysis. Parameters (α, β, K k , K P , s) that define the kinetics of U, V were fixed to those used earlier [70]. Two to four independent simulations from different random seeds were executed for a given parameter set to obtain average feature scores F. Because grid-search in 8 dimensional parameter space was unattainable due to heavy computational load, we adopted parameter search around manually selected reference parameter sets. First, based on preliminary simulations performed with various parameters, a single parameter set R1 was chosen which gave rise to polarized cell morphology; i.e. elongated shape in the moving direction. Around the reference set R1, the following parameters were varied: a w and W tot which appeared to affect the elongation of cell in the moving direction,χ U and D W related to the split of the leading edge, and γ that seems to affect the appearance of new pseudopods. We performed a grid search around R1 in 2-dimension (a W , γ) for χ U = 0 and χ U = 50 (S5A Fig left panels) while fixing other 6 parameters. Similarly, a grid search was performed in (W tot , D W ) for χ U = 0 and 50 ( S5A Fig left panels). Second, we chose another parameter set R2 which gave rise to relatively round shape and random cell displacement (i.e. near HL60-like behavior). Around R2, we varied (a W , W tot ), (γ, ρ), and (k W1 , μ) at both χ U = 0 and 50 (S5A Fig right panels)

Pseudopod analysis
For both experiments and simulation data, events of pseudopod formation were detected manually by eye as either 'de novo', 'Y-split' or 'one-way-split' based on [7]. New pseudopods that appeared in non-protruding regions (Fig 4E left) were counted as "de novo formation". Protrusions that appeared close to the leading edge were marked as 'Y-split' because of the resulting cell shape (Fig 4E right). When new protrusive events occurred on top of the leading edge so as to offset and steer the extending pseudopod, it was counted as 'one-way-split' (Fig 4E  middle) [7]. For the subset of the data, the analysis was performed three independent times to confirm reproducibility. To obtain angular distribution, direction of pseudopod extension 10 seconds from the onset was measured relative to the centroid movement. Aggregation-stage Dicyostelium (N = 90 from 4 timeseries), and the closest simulations (N = 44 from 2 timeseries). HL-60 (N = 68 from 5 timeseries), and the closest simulations (N = 77 from 3 timeseries).
Supporting information S1 Text. The online supporting information includes. 1. Details of our mathematical model.