Random walk and cell morphology dynamics in Naegleria gruberi

1 Introduction

Combinatorial use of persistent motion and reorientation is a common feature found in cell movement. Be it bacterial swimming or amoeboid crawling, persistent movement allows cells to gain most distance in one preferred direction so as to facilitate efficient escape from hazards or conversely attraction to nutrients. Reorientation on the other hand is not only required to adjust direction of persistent movement but also to facilitate cells to randomly explore and survey uncertain extracellular environments (Viswanathan et al., 1999; Bartumeus et al., 2002). In E. coli bacteria, the cell movement consists of a period of straight run interrupted by a stall or “tumble” where flagellar rotation reverses and cells reorient in random directions. The frequency of tumbling is regulated through a chemosensory system so as to provide orientation bias towards an attractant or away from a repellent. The exact nature of such motility pattern determines how well E. coli cells disperse (Taktikos et al., 2013). In the amoeboid movement, pseudopodal protrusions enriched in branched F-actin networks (Pollard, 2007) are formed transiently and can guide cells in different orientations. In addition, a confined region of the plasma membrane needs to retract in order to realize net displacement. In many cell types, cortical F-actin that is crosslinked with myosin II is enriched in such contractile membrane regions (Chi et al., 2014). Persistent movement arises when a cell has mono-polarity meaning that it has a single dominating leading edge and a retracting trailing end. The occurrence and location of these organizational events along the plasma membrane determine the sequential appearance of plasma membrane protrusions and rear retractions, ultimately influencing the direction, speed, and duration of cell movements.

Quantitative time-series analyses of cell displacement and cell shape change are important for explicit characterization of random cell motion. In many cases, cell displacement can be approximated as a particle obeying persistent random walk. Phenomenologically, the simplest form of differential equation that describes such stochastic dynamics is the Langevin equation (Dunn and Brown, 1987; Selmeczi et al., 2005; Selmeczi et al., 2008)

$\frac{d\overrightarrow{v}\left(t\right)}{dt}=-\beta \overrightarrow{v}\left(t\right)+\sigma {\overrightarrow{\xi }}_t(1)$

where $\overrightarrow{v}$ is the velocity vector, β is the decay rate, σ is the noise strength, and ${\overrightarrow{\xi }}_t$ is 2D white Gaussian noise. Random walk of E. coli can be approximated by Brownian motion having a short-term memory. In eukaryotic crawling, cell trajectories of fibroblast cells (Gail and Boone, 1970) and endothelial cells (Stokes et al., 1991) are also consistent with this simple persistent random walk model. In many other cell types, random walk includes memory that depends on the velocity and orientation which can be described by modifications to the above model (Takagi et al., 2008; Li et al., 2011). There are also random walk statistics called Lévy-walk with step lengths that follows a long tail (power-law) distribution (Viswanathan et al., 1999). There, the Mean Square Displacement (MSD) essentially diverges and the trajectories are characterized by self-similarity of the step lengths (Reynolds, 2010; Reynolds and Ouellette, 2016). Because Lévy-walk has very small probability of revisiting the same location, it is thought to arise in systems such as bird foraging that require an efficient search strategy. At the cellular-level, effector T-cells (Harris et al., 2012) and swarming bacteria have been reported to exhibit Lévy-flight like statistics (Ariel et al., 2015; Huo et al., 2021).

To date, quantitative understanding of random walk behavior of amoeboid cells is limited to data from a handful of cell-types; these are mostly timelapse microscopy images of cultured metazoan cells and amoebozoa Dictyostelium. From microbial ethology and evolutionary biology perspectives, however, we should note that amoeboid movement is found not only in animals, fungi and amoebozoans (Prostak et al., 2021) but also in largely uncharacterized species in the excavate taxon Heterolobosea namely Naegleria spp. and the slime mold acrasids (Brown et al., 2012). The ancestors of the opisthokont lineage and Naegleria diverged more than a billion years ago (Parfrey et al., 2011). Knowing the details of motility characteristics in Naegleria should help us understand the common design of the motility trait that is either deeply conserved across taxa or acquired independently by strong selective advantages.

Among members of genus Naegleria, non-pathogenic Naegleria gruberi (hereafter refer to as N. gruberi) is the better characterized species whose genome has been sequenced (Fritz-Laylin et al., 2010). In its amoebic phase, N. gruberi grows and divides by feeding on bacteria through phagocytosis (Fulton, 1970). Under low electrolyte conditions, it quickly shifts to the non-feeding flagellated state by rapid de novo synthesis of microtubules (Walsh, 2007). In the amoebic state, the overall cell cortex is enriched in F-actin with marked accumulation around membrane ruffles (Velle and Fritz-Laylin, 2020). An early work using reflection interference microscopy has revealed that N. gruberi adhere and form discrete dot-like contacts to non-treated glass surfaces and migrate (Preston and King, 1978). These so-called “focal contacts” leave behind footprints of membrane residues on the glass substrate as the cells crawl away (Preston and King, 1978). With the advent of genomics and molecular cell biology, it has become clear that N. gruberi possess the essential side-branching nucleator of F-actin—the Arp2/3 complex and its activators WASP and SCAR (Fritz-Laylin et al., 2017; Velle and Fritz-Laylin, 2020; Prostak et al., 2021). Inhibition of Formin reduces directional persistence, and inhibition of the Arp2/3 complex reduces the cell speed (Velle and Fritz-Laylin, 2020). N. gruberi also has Myosin II (Sebé-Pedrós et al., 2014) and a potential orthologue of Integrin beta (Morales et al., 2022), although whether they exist in other groups in Excavata is unclear (Velle and Fritz-Laylin, 2019).

While the above works indicate likely similarity of actin-dependent processes involved in cell crawling in an evolutionary distant eukaryote, quantitative characterization of the cell-level motility pattern is so far lacking. Do N. gruberi cells exhibit persistent random walk behavior? What is the characteristic time scale of persistence and reorientation if any? How similar are their movements compared to the well-studied systems such as Dictyostelium and immune cells? In this report, we performed quantitative analysis on cell movements and shape change of N. gruberi. Our analysis demonstrates that N. gruberi cells exhibit persistent random walk driven by a large morphology change that involves appearance, splitting and annihilation of uniquely complex pseudopodium protrusions.

2 Results

2.1 An overview of cell movement and cell morphology

To quantitate the movement of N. gruberi on a two-dimensional flat surface, we performed phase contrast time-lapse microscopy. A non-coated glass coverslip was employed as a cell substrate throughout this study. Figure 1A shows representative phase contrast images of N. gruberi in liquid growth media (Materials and Methods). The cells under our culture condition exhibited one or more hyaline protrusions that appeared dark in phase contrast images (Figure 1A arrows). In the example shown, protrusions extended along the glass surface for 15–50 s and the one that became dominant (i.e. the leading edge) extended in the direction of the overall cell movement (Figure 1A, 0 s). Marked cytoplasmic streaming from the center of the cell towards these extensions was observed (Supplementary Movie S1). A new protrusion appeared and extended first in the lateral direction (Figure 1A 20 s, 60 s arrow) and steered towards the front. It was then bent sideway before being retracted (Figure 1A 40 s, 80 s). Duration of the pseudopod extension/retraction cycle varied between 15–50 s (Supplementary Figure S1; Supplementary Movie S2). Concomitant reversal of cytoplasmic streaming was observed during retraction of pseudopods. A small bud-like bulge at the trailing end of a cell which we shall refer to as “uroid” appeared as a residue of a retracted pseudopod that was retained for an extended period of time (Figure 1A white circle). The uroid contained thin filopodia-like projections as described earlier (Preston and King, 1978).

FIGURE 1

FIGURE 1. An overview of N. gruberi movement. (A) Representative phase contrast images from a time-series of a migrating N. gruberi cell. Arrows: protruding edges. Circles: a bud-like rear structure (“uroid”). (B) First 360 s of randomly selected centroid trajectories. 4 trajectories are separately shown for visibility. Scale bar: 100 μm.

Under our culture condition, the cells appeared to re-orient in random directions at irregular timing. We performed long time cell tracking by employing an automated stage that was programmed to track target cells (see Methods). Figure 1B shows representative cell trajectories obtained from the automated tracking. The trajectories consisted of a period of straight movement that lasted for about 30–200 s and a time period of relative low displacement and re-orientation (Figure 1B). The movement is thus, at surface, akin to the run-and-tumble behavior of E. coli. There was a close link between the run/re-orientation dynamics with the cell shape. During a straight run, cells took a fan-like shape (Figure 2A; Supplementary Movie S3). The tail remained narrow while the front was occupied by a broad lamellipodia that expanded then split into branches of pseudopods (Figure 2A, 0–16 s). These bifurcating protrusions often fused to restore a large lamellar extension (Figure 2A, 20 s). On the other hand, cells re-oriented when the bifurcated protrusions remained separate. In most cases, the uroid persisted during front splitting and thus the cells took a Y- or trident shape (Figure 2B; Supplementary Movie S4). There were also cases where the uroid disappeared in Y-shaped cells (Figure 2C; Supplementary Movie S5). The two fronts expanded in the opposing directions and gave rise to a transient “dumbbell-like” bipolar morphology (Figure 2C, 70 s). After 10 s, one end shrunk and became the uroid while the other end became the next front (Figure 2C, 80 s). There was little centroid displacement during this period which lasted for about 40 s.

FIGURE 2

FIGURE 2. Protrusion dynamics and the cell shape change. (A) A fan-shaped cell with front splitting during a persistent run. (B) Front splitting followed by reorientation (curvature kymograph for the sequence is shown in Figure 5A). (C) Dumbbell shape arise after front splitting and disappearance of the uroid (curvature kymograph for the sequence is shown in Supplementary Figure S4E). Arrows: orientation of centroid movement. Inverted triangles: propagating curvature waves.

2.2 Random walk statistics

To characterize the random-walk statistics, we quantified the mean square displacement (MSD) and the instantaneous speed defined by the centroid displacement in 1 s time interval from trajectories of N = 10 cells (Figure 3). Even with the help of automated stage tracking, fast movement of N. gruberi made it difficult to track cells for long time duration before they come close to the edge of the plate or collided with one another. Thus, to obtain MSD, single trajectories were each divided into sub-trajectories of 100–3,600 s time-window and treated as independent data samples (Figure 3A). The slope of the MSD from the individual trajectories was 1.5–2.0, where the mean and standard deviation are 1.77 and 0.08 (Figure 3B; Supplementary Figure S2A). The time-dependency of the MSD indicates that the random walk of N. gruberi falls somewhere between pure Brownian (exponent of 1) and ballistic (constant velocity) motion (exponent of 2). Figure 3C and Supplementary Figures S2B–F show the distribution of the instantaneous velocity. The distribution followed 2-dimensional Gaussian with zero-mean and standard deviation of 51 μm/min (0.86 μm/s) (Figure 3C, Supplementary Figures S2B–F). This feature is distinct from Dictyostelium random motility which is non-Gaussian (Takagi et al., 2008). The median of the absolute speed was 60 μm/min (1.0 μm/s) which is close to the average speed reported in earlier literatures (King et al., 1981; Thong and Ferrante, 1986).

FIGURE 3

FIGURE 3. Statistics of the centroid movement. (A) The trajectories used for detailed analysis (N = 35). The origin is set at the initial position. Scale bar: 500 μm. (B) The ensemble averaged MSD. Red solid curve line: exponent 1.8, which is the result of fitting the ensemble averaged MSD. Shaded region: standard deviation. (C) Probability distribution of $\left|\widehat{\overrightarrow{v}}\right|$. Red solid curve: a fitting curve from a 2-dimensional isotropic Gaussian distribution with a standard deviation of 0.86 μm/s. (D) The ensemble averaged VAC. Red solid curve: the sum of two exponential functions (Eq. 2). (E) Acceleration parallel (blue) and orthogonal (green) to the velocity. The binned average (circle) and the standard deviation (star). (F) Persistence of displacement orientation in a unit time step. Cosine of the angle change θ in the velocity in 1 s interval. The binned average (circles) and the standard deviation (error bars).

The temporal auto-correlation of the centroid speed (velocity auto-correlation; VAC) (Supplementary Figure S2G) shows, on average, that there are two characteristic decay times that cross over at around 10 s (Figure 3D). By assuming that VAC follows the sum of two exponential function (Selmeczi et al., 2005) with velocity $\overrightarrow{v}$:

${〈\overrightarrow{v}\left(t+\tau \right)\bullet \overrightarrow{v}\left(t\right)〉}_t={\mathrm{\Phi }}_1{e}^{-\frac{\tau }{T_1}}+{\mathrm{\Phi }}_2{e}^{-\frac{\tau }{T_2}}(2)$

we obtained decay time T₁ and T₂ of approximately 6 s and 90 s, respectively, where the weight ${\mathrm{\Phi }}_1$ and ${\mathrm{\Phi }}_2$ are 0.36 μm²/s² and 0.87 μm²/s² (Figure 3D Red curve). Based on the Bayesian information criterion (BIC) (Schwarz, 1978), two exponential functions gave the lowest value compared to one or three (Supplementary Figure S2H). When the length of time sequence chosen for the analysis was doubled from 50 s to 100 s, deviation of the parameter values was within an order of magnitude (Supplementary Figure S2H magenta curve). Decay time T₂ of approximately 90 s was also evident from the time derivative of VAC (Supplementary Figure S2I).

To check for orientational preference in the memory, we plotted the relationship between velocity and acceleration (change in $\overrightarrow{v}$ in 1 s interval) (Supplementary Figures S2J, K). The mean acceleration orthogonal ($a_{\perp }$) to the velocity was near zero regardless of |v| (Figure 3E; green circle) with non-zero variance (Figure 3E green stars) which suggests that the orientation of N. gruberi has no apparent left-right asymmetry. On the other hand, the mean acceleration parallel ($a_{\parallel }$) to the velocity was near zero at small velocity then decreased towards the negative at large velocities (Figure 3E). The standard deviation (Figure 3E; blue stars) increased somewhat at high |v|, however rarity of these fast step events prevented us from obtaining reliable averages. These features of acceleration are similar to those reported for Dictyostelium (Takagi et al., 2008). The negative acceleration parallel to the migration direction at high |v| implies that the cells do not maintain high |v| during re-orientation. Accordingly, when we plot reorientation angle θ as a function of |v| (Figure 3F) we see that most of re-orientation occurs below |v| = 1 μm/sec. Above 1 μm/s the cells are moving in a straight line; i.e. cosθ = 1.

2.3 Generalized Langevin equation

To gain further insights on the specifics of the random walk statistics, it is instructive to compare the data with the behavior of simple idealized equations. The velocity auto-correlation that follows the sum of two exponential indicates that random walk dynamics cannot be captured simply by the Ornstein-Uhlenbeck process (Eq. 1) which only has a single exponent (Dunn and Brown, 1987). A straight-forward and minimal extension to Eq. 1 is to include additional memory with the decay rate γ as an integral in the form of generalized Langevin-equation (Selmeczi et al., 2005)

$\frac{d\overrightarrow{v}\left(t\right)}{dt}=-\beta \overrightarrow{v}\left(t\right)+{\alpha }^2{\int }_{-\infty }^t{e}^{-\gamma \left(t-t^{\prime }\right)}\overrightarrow{v}\left(t^{\prime }\right)d{t}^{\prime }+\sigma {\overrightarrow{\xi }}_t(3)$

Here, α is the strength of memory effect, and ${\overrightarrow{\xi }}_t$ is a normalized Gaussian white noise that satisfies $〈{\overrightarrow{\xi }}_t〉=0,〈{\overrightarrow{\xi }}_t{\overrightarrow{\xi }}_{t^{\prime }}^{\mathrm{T}}〉=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\delta \left(t-t^{\prime }\right)$,$〈〉$ is an ensemble average and $\mathrm{\delta }\left(t\right)$ is the delta function, σ is the noise strength (Selmeczi et al., 2005). By introducing

$\overrightarrow{V}\left(t\right)=\alpha {\int }_{-\infty }^t{e}^{-\gamma \left(t-t^{\prime }\right)}\overrightarrow{v}\left(t^{\prime }\right)d{t}^{\prime }(4)$

the equation of motion becomes

$\frac{d\overrightarrow{v}\left(t\right)}{dt}=-\beta \overrightarrow{v}\left(t\right)+\alpha \overrightarrow{V}\left(t\right)+\sigma {\overrightarrow{\xi }}_t(\mathrm{5a})$

$\frac{d\overrightarrow{V}\left(t\right)}{dt}=\alpha \overrightarrow{v}\left(t\right)-\gamma \overrightarrow{V}\left(t\right)(\mathrm{5b})$

Based on the values of $T_1,T_2,{\mathrm{\Phi }}_1,{\mathrm{\Phi }}_2$ obtained above, we calculated the parameter values of the generalized Langevin equation (Eqs 5a, b) from the analytically obtained VAC at the steady state (see Eq. 29).

Trajectories, the MSD and the VAC were obtained by numerically calculating Eqs 5a, b with the parameters obtained above (Table 1). The individual trajectories consist of combination of persistent movement and turns (Figure 4A). The slope of MSD had mean and standard deviation of 1.80 ± 0.06, which matched well with the experimental data (Figure 4B). The distribution of |v| showed a single peak that was slightly smaller compared to the experimental data (Figure 4C). The median was 56 μm/min (0.94 μm/s) in the simulation, which matched well with 60 μm/min in the experiment. The velocity autocorrelation consists of two slopes that crossed over at around 10 s (Figure 4D red), which was similar to the crossover in the experimental data (Figure 4D black). Velocity dependence of acceleration also matched well with the experimental data (Figure 4E). On the other hand, the range of cell speed at which turning occurred in the simulations was somewhat broader (0–1.2 μm/s) compared to the real cell (0–0.8 μm/s) (Figure 4F). While the MSD and the VAC characteristics were well captured by the memory effect described in Eq. 3, deviation from the model became evident when comparing autocorrelation separately for the centroid movement (absolute velocity $\left|\overrightarrow{v}\right|$) and the orientation ($\overrightarrow{v}/\left|\overrightarrow{v}\right|$) (Supplementary Figure S3). In the experimental data, it is only the autocorrelation of the orientation $\overrightarrow{v}/\left|\overrightarrow{v}\right|$ not $\left|\overrightarrow{v}\right|$ that showed two decay times (Supplementary Figure S3A, B). In the generalized Langevin-equation, the velocity and the orientation share the same time scales, and thus the autocorrelation of both the orientation $\overrightarrow{v}/\left|\overrightarrow{v}\right|$ and $\left|\overrightarrow{v}\right|$ decayed with the two exponents (Supplementary Figures S3C, D).

TABLE 1

TABLE 1. Parameters for the generalized Langevin equation. The experimental data were fitted with the analytical VAC (Eq. 29).

FIGURE 4

FIGURE 4. Statistics of the persistent random walker trajectories. (A) Simulated trajectories. (B) MSD; the ensemble average (circle) and the standard deviation (shade) (red). Experimental data (black) are duplicated from Figure 3B for comparison. (C,D) Probability distribution of $\left|\widehat{\overrightarrow{v}}\right|$ (C) and VAC (D) (red). Experimental data from Figures 3C, D (black) are duplicated for comparison. (E) Acceleration, parallel (red) and orthogonal (black) to the velocity. The average (cross mark) and the standard deviation (triangle). Blue and green markers show the experimental data in Figure 3E. (F) The average of cos θ (red circle) and the standard deviation (red error bar). Experimental data from Figure 3F (black circle and error bar) are duplicated for comparison.

2.4 Cell shape dynamics

Rather than pursuing extensions of the particle-based formalism such as those that treat the two timescales separately (Li et al., 2008; Takagi et al., 2008), we sought to more directly characterize cell reorientation by analyzing the cell shape dynamics. Based on binarized cell mask images and a boundary tracking algorithm (Nakajima et al., 2016; see also Supplementary Figure S4A), 500 points along the edge of cell masks were tracked in the laboratory frame for the local curvature and the normal velocity (Figure 5A, B; see also Supplementary Figures S4B, D, F for another sample). A protruding edge can be seen as a positive local-maximum in the curvature (Figure 5A yellow regions). The advancing front of a cell can be discerned by its positive velocity (Figure 5B, yellow regions), and the trailing uroid by the negative velocity (Figure 5B, blue regions). At the cell front, a new protrusion frequently appeared to split off from a pre-existing protrusion (Figures 5A, B white arrows). These appeared in the kymograph as branching positive curvature regions that propagated rearward until they were annihilated at or near the uroid (Figures 5A, B black arrows). The sequence of curvature wave dynamics represents well the shape dynamics as seen in the snapshots (Figures 5C, D orange arrows; see also Figure 2B white arrows for a protrusion from split to annihilate).

FIGURE 5

FIGURE 5. Cell boundary analysis. (A,B) The curvature (A) and the normal velocity (B) of the cell boundary taken from a representative cell exhibiting random walk. White arrows: splitting. Black arrows: pair annihilation. (C,D) Magnified view of a subsection in (A,B). Orange arrows indicate protrusions that split (C) or annihilated (D). (E,F) The curvature (E) and the normal velocity (F) of the boundary taken from a cell with high persistence. (G) Snapshots of the cell analyzed in (E,F). The white arrow: the direction of the centroid movement. The inverted triangles mark the uroid.

The curvature wave dynamics are surprisingly similar to those obtained for Dictyostelium and neutrophil-like HL60 cells (Driscoll et al., 2012; Imoto et al., 2021) with a noticeable difference that splitting was more frequent and thus numerous. The other difference compared to Dictyostelium and HL60 cells is the occasional and transient appearance of dumbbell-like shape (Supplementary Figures S4C, E, G; Supplementary Movie S5). When it appears, the centroid velocity orientation showed discontinuous change (Supplementary Figure S4C, black arrow). In the kymograph representation, a dumbbell-like cell shape appears as two or three stable curvature waves (Supplementary Figure S4E, black arrow). Most positions had zero velocity (Supplementary Figure S4G, black arrow), indicating stalling of cell shape change. These observations indicate that as the dumbbell shape appeared, the cell stopped and randomized its orientation. There were also rare cases where the cell maintained mono-polarity for an extended period of time (Figures 5E–G; Supplementary Movie S6; see Supplementary Figure S5 for additional samples). There, new curvature waves emerged frequently and traveled fast before disappearing at the tail (Figure 5E). The position where curvature waves appeared always showed positive velocity, while the positions where curvature waves disappeared showed negative velocity (Figures 5E, F). These patterns in the kymograph correspond well with the observation of fast curvature waves that propagate from the advancing cell front and disappear at the uroid (Figure 5G).

A further analysis showed a close relationship between the curvature wave and the centroid movements. The protruding and the retracting membrane regions can be identified as positive curvature regions with positive (Figure 6A, white dots) or negative (Figure 6A black dots) velocity respectively. The orientation of the normal vector at the protruding region showed high correlation with the direction of centroid velocity (Figure 6B blue). The retracting regions oriented in the opposite directions which appeared somewhat broader in distribution (Figure 6B orange). To further analyze the dynamics of the curvature wave, high curvature regions (Figure 6C white) at each time frame were assigned as individual protrusions (Figure 6C green). While there were multiple protrusions in the protruding region, a dominant leading edge can be detected from identifying a single protrusion whose normal vector angle was the closest to that of the centroid velocity (Figure 6C magenta). Once a curvature wave became the leading edge, it remained so for about 2.8 s as measured from its average lifetime (Figure 6D). Another interesting feature of the membrane extensions is that they gave birth to secondary pseudopods or were steered to other directions. The typical angular velocity associated with this dynamic was 0.1 rad/s (Figure 6E). Together with the two timescales of decay (Supplementary Figure S3B), these behaviors indicate that the centroid velocity angle by itself follows 1D persistent random walk. From experimentally obtained parameters of the leading edge lifetime (2.8 s) and the angular velocity 0.1 rad/s, the 1D model (see Materials and Methods, Cell Boundary Analysis section) yielded decay time of 142 s on average which matched well with the experimental data (Figure 6F).

FIGURE 6

FIGURE 6. Relation between the membrane protrusions and the centroid velocity angle. (A) The protruding front (white) and the retracting rear (black) detected from the velocity kymograph are overlaid on top of the curvature kymograph (see Methods). (B) The angular histogram of the protruding front (blue) and the retracting rear (orange) relative to the cell orientation as determined by the centroid velocity. (C) The position of protrusive regions (“curvature waves;” green). The region that co-extended most closely in the direction of the cell centroid motion (“leading protrusion;” magenta). The binarized mask of the protrusion region (white) obtained from the curvature kymograph is shown in the background. (D) Duration time histogram of the leading protrusion (magenta in (C)). (E) Histogram of the angular velocity in the protrusion orientation [the vector normal to the cell contour at positions indicated in green in (C)]. Solid lines are exponential fit to the data (D,E). (F) Estimated VAC decay time $T^{est}$ for the representative data.

2.5 Fourier-based morphology space analysis

To obtain a quantitative morphometry, we chose by eye 21 representative mask images each for the 3 shapes; fan-shape, split and dumbbell (Supplementary Figure S6A) and calculated the Fourier power spectrum of the cell edge coordinates and their principal components were calculated (see Materials and Methods). We found that the first two principal components were sufficient to obtain well separated clusters that represented the shape class (Figure 7A). All cell masks analyzed were distributed within a confined domain in the PC1_fourier -PC2_fourier space (Supplementary Figure S6B). The fan-shaped data were located at a low PC1_fourier and high PC2_fourier region (Figure 7A circles). The split-shape were found in the low PC1_fourier—low PC2_fourier region (Figure 7A asterisks). The dumbbell-shape was located at high PC1_fourier and high PC2_fourier (Figure 7A triangles). To see what shape features the principal components represented, we reverse calculated an artificial form by obtaining Fourier spectrum as a product of synthetic principal component vector to the eigen vector matrix composed of the basis of Fourier spectrum (see Methods). In brief, PC1_fourier indicated the aspect ratio i.e., elongation, PC2_fourier the head width, PC3_fourier the rear steepness (Supplementary Figure S6C). Here, the main contribution to PC1 were from the wave number 1 and −1 with coefficients of 0.68 and 0.73. For PC2, the contribution from wave number 1, −1, 2, and 3 was 0.62, −0.59, −0.49, −0.12, respectively. Contribution from other modes was small with coefficients less than 0.03.

FIGURE 7

FIGURE 7. Fourier analysis of the cell contour. (A) Principal component space (PC1_fourier, PC2_fourier) obtained from 63 manually selected binarized snapshots (left panel). Representative cell masks (right panels). (B) Time series in PC1_fourier-PC2_fourier space from 3 representative timelapse sequences (colors). The time spent in the negative PC1_fourier region per total trajectory time was 676 s/814 s (blue), 287 s/382 s (orange), and 459 s/578 s (green). (C) Time evolution of PC1_fourier and PC2_fourier of 10 s around a large turn that involves transition to the dumbbell shape (3 representative events; colors). Black arrow indicates the direction of time evolution. (D) Time evolution of PC1_fourier and PC2_fourier during persistent migration. Colors indicate different time series [duration: 269 s (blue), 1,039 s (orange), or 3,600 s (green)]. (E) Autocorrelation of C_n (n = $\pm$ 2, $\pm$ 3, $\pm$ 4). The decay rate: 12.2 (C₂), 18.0 (C_-2), 8.9 (C₃), 7.6 (C_-3), 10.6 (C₄), and 4.4 (C_-4) seconds. (F) Distribution of angles of centroid velocity and $C_{34}$. (G,H) The x- (G) and y-components (H) of centroid velocity plotted against real (G) and imaginary (H) parts of C₃₄. Red circles and blue bars indicate average and standard deviation of centroid velocity binned with the value of C₃₄. Orange lines indicate the result of fitting with linear proportionality.

How the cell shape changed during turning can be analyzed by tracking the time sequence in the PC1_fourier—PC2_fourier space. Figure 7B shows three independent samples of re-orienting cells. Here, cells were mainly located in the negative PC1_fourier region with occasional visits to the positive PC1_fourier region. This is consistent with the above observation that cells took fan- or branched-shape (negative PC1_fourier) in addition to rare occurrence of dumbbell-shape (positive PC1_fourier). Figure 7C shows three independent samples of the dumbbell-shape forming cells. The counter-clockwise circular trajectories in the PC1_fourier—PC2_fourier space signify a transition from the fan-shape to splitting then to the dumbbell-shape. On the other hand, Figure 7D shows three independent trajectories that remained in the negative PC1 region for extended period of time. These cells at least during the time window of observation fluctuated between the fan-shape and the bifurcating fingers. There was no clear relationship between the morphometry state (PC1_fourier, PC2_fourier) and the cell speed (Supplementary Figures S7A, B). There was, however, negative correlation between the centroid speed and the rate of state transition $d\left\{\text{PC}{1}_{\text{fourier}}\right\}/dt$ but not with $d\left\{\text{PC}{2}_{\text{fourier}}\right\}/dt$ (Supplementary Figures S7C, D). As the former relation was seen in the negative direction $d\left\{\text{PC}{1}_{\text{fourier}}\right\}/dt$ <0, it signifies that cells accelerate when recovering from dumbbell-shape.

Besides the rate of state transition in the principal components, there should be a direct relationship between the Fourier components C_n themselves and the centroid movement. Autocorrelation analysis showed that the decay rates for C_-3, C₃, and C₄ were 7.6, 8.9, and 10.6 s, respectively (Figure 7E) and thus matched most closely to the short decay time of VAC. As for the centroid velocity itself, according to the deformation tensor-based theory of cell movement (Ohta et al., 2016), it should be proportional to $C_{nm}\equiv {\dot{C}}_{-n}{C}_m-C_{-n}{\dot{C}}_m$ where −n + m = 1. More specifically, C_nm is a complex number whose absolute value |C_nm| and the angle arg (C_nm) are expected to be proportional to the speed and the velocity angle of the centroid respectively. In NIH3T3 cells, it has been shown that velocity is proportional to the elongation C_-2 and triangular C₃ modes of deformation multiplied by their time derivatives; i.e. $C_{23}={\dot{C}}_{-2}{C}_3-C_{-2}{\dot{C}}_3$ (Ebata et al., 2018). However, in N. gruberi, we found little correlation between C₂₃ and the centroid velocity (Supplementary Figures S7E–G). Instead, we found that it was $C_{34}={\dot{C}}_{-3}{C}_4-C_{-3}{\dot{C}}_4$ that correlated highly with the centroid velocity angle (Figure 7F) and x- and the y-component of the centroid velocity (Figures 7G, H). The difference between Naegleria and NIH3T3 may be attributed to the fact that Naegleria has many pseudopods that are complex in shape as analyzed below.

2.6 Deep learning-based morphology analysis

To further investigate the cell shape characteristics, we employed a convolutional neural network that was previously trained to classify cell shapes based on similarity to Dictyostelium-like, HL60-like, or fish keratocyte-like shapes (Imoto et al., 2021). While the method is not suited to track shape change over time due to discrete change in the morphometry space that is sometimes introduced by uncertainty in the cell orientation during mask alignment, it has an advantage of providing an objective morphometry that is independent of known feature basis. On average, Naegleria was classified as Dictyostelium-like (high PC1_cnn, low PC2_cnn) (Figure 8A). This was natural as it has been shown to pick up branching shapes that are elongated overall in the migrating direction (Imoto et al., 2021). We noticed substantial variability, however, in the individual cell shapes (Figure 8B; black) that exceeded those normally observed in Dictyostelium (Figure 8B; green). Shapes that deviated in the PC1_cnn direction were mapped to dumbbell-like domain in the Fourier descriptor-based morphometry (Figure 8C; orange). Those that deviated towards low PC1_cnn were mapped to the domain that showed numerous pseudopods (Figure 8D; orange). Datapoints that fell at or near the HL60-like domain (low PC1_cnn, low PC2_cnn) were mostly fan-like (Figure 8C; blue; Figure 8D; blue) and their occurrence per timeseries showed positive correlation with the MSD exponent (Figure 8E). This is consistent with the notion that more mono-polarized the cells are, the more ballistic the cell trajectories become.

FIGURE 8

FIGURE 8. Shape analysis by a CNN-based classifier. (A,B) Time-average (A) or snapshot (B) of PC1_CNN and PC2_CNN from N. gruberi images (black) were superimposed on the PC1_CNN -PC2_CNN phase space of D. discoideum (green), HL60 (red), and fish keratocyte (yellow). (C) Snapshot data in (B; black) was classified into split (black), fan (blue) or dumbbell (orange) based on PC1_fourier and PC2_fourier. (D) PC1_fourier and PC2_fourier of HL60-like (blue), cells with PC1_CNN lower than −1.6×10⁶ (orange), or the other cells (black) classified with CNN. (E) Ratio of frames whose shape was classified as HL60 in deep leaning-based classification.

3 Discussion

In this report, we analyzed movements of N. gruberi cells by quantifying their speed, directionality, and shape change. The locomotive speed of N. gruberi cells was around 60 μm/min, which is similar to that reported in early literatures (King et al., 1981; Thong and Ferrante, 1986). It is substantially larger in magnitude compared to that of fibroblast ∼0.4–1.0 um/min (Welf et al., 2012; Passucci et al., 2019), and even larger compared to fast migrating cells such as vegetative Dictyostelium 5 μm/min (Li et al., 2008), and neutrophils 17 μm/min (Hartman et al., 1994). Despite the large speed difference, we found common features between N. gruberi and other cell types whose random motility have previously been characterized. The exponent of MSD was approximately 1.8 meaning that the random walk is non-Fickian and non-ballistic at least at surface. Stronger deceleration at higher velocity implies non-ballistic movement, where the non-memory term, i.e., fluctuating components plays a dominant role in determining the next move. Similar exponent is known in MDCK cells (Dieterich et al., 2008), A549 cancer cells (Kwon et al., 2019) hematopoietic progenitor cells (Partridge et al., 2022), and T cells (Jerison and Quake, 2020). Of particular note is that the time-scale where such exponent was observed for N. gruberi was about 10–100 s which is within the order of magnitude required for a cell to move one cell-body length. This seems also to be the case for MDCK cells where the exponent of 1.8 was observed at much longer time-scale of 4–20 min with corresponding length scale 4 μm–20 μm. All in all, our data combined with the observations above from earlier literatures suggest that the time scale at which cells move in a straight line is the major determinant of cells’ displacement.

The other common feature found in this study was the presence of two characteristic decay time in the VAC (Selmeczi et al., 2008). For N. gruberi, these were T₁ = 6 and T₂ = 90 s, which are in the same order of magnitude as that of Dictyostelium in the vegetative (T₁ = 5.2 and T₂ = 228 s) and the starved (T₁ = 11 and T₂ = 108 s) states (Li et al., 2008). Although equivalent measurements have not been documented for neutrophils, their cell shape changes had typical time scale of 8 s (Hartman et al., 1994) and the persistence time during chemotaxis was 103 s (Itakura et al., 2013; Haastert, 2021). From the MSD measurement, the persistent time of Dictyostelium and neutrophil-like HL60 were 151 s and 278 s, respectively (Imoto et al., 2021). Interestingly, VAC of Human keratinocyte-like cells (HeCaT) whose speed was much slower (0.18 μm/min) could also be fit with the sum of two exponential functions [T₁ = 76 s and T₂ = 860 s; (Selmeczi et al., 2005)]. The characteristic time scale of around 10 s was attributed to the time scale of actin polymerization in the protruding pseudopodia (Haastert, 2021). However, the pseudopod lifetime in N. gruberi was rather long; about 15–50 s. The discrepancy may be attributed to the sister pseudopods that formed from the main pseudopods which were not analyzed in our manual tracking. In some cases, the pseudopod itself also appeared to bend in one direction. In support of this notion, the autocorrelation of C_-3 and C₄ had decay time of 7–10 s which matched well with the first decay time of VAC.

On the other hand, the second decay time of VAC (90 s; Figure 3D) was close to the timescale of directional persistence i.e. “run” phase estimated from the curvature wave dynamics (142 s; Figure 6F). As for the average cell speed, we found strong correlation between the centroid velocity with the coupling of deformation modes C_-3 and C₄, instead of C_-2 and C₃. This suggests that the orientation of N. gruberi cells depends not on the primary membrane protrusions but on their sister sub-structures. A further pseudopod-level analysis at finer time-scale is required to clarify the relation between the deformation modes and the branching pseudopods. The rare cells with high persistency did not take high PC1_fourier value (Figure 7D) which was opposite of Dictyostelium (Tweedy et al., 2013). This likely stems from the fact that, in N. gruberi, the elongated form was usually dumbbell-shaped which occurred when the cells stalled and reoriented.

The splitting pseudopod may entail a mechanism similar to those found in amoebozoan and metazoan cells where dendritic actin meshworks are regulated by excitable and oscillatory dynamics (Huang et al., 2013). The presence of local inhibitor of pseudopod formation in neutrophils and Dictyostelium (Xu et al., 2017) and potential lack of such in Naegleria may underlie the difference in the number of pseudopods. Alternatively, there may be local reduction in the actin cortex that are stochastic in nature. Although protrusions observed under our culture conditions did not appear as blebs, marked flow of cytosol towards the membrane observed during extension of a protrusion suggests local pressure release. Protruding form triggered by the pressure difference at a fluid-fluid interface is known as viscous fingering. The movement speed of N. gruberi was 5 times as large as that of neutrophils and Dictysotelium, but closer to that known for fragments of Physarum which also exhibit marked cytoplasmic streaming (Rieu et al., 2015) and persistent random walk (Rodiek and Hauser, 2015). Such high velocity and potential interface instability may underlie the observed branching of pseudopods. Another unique shape feature was the dumbbell-like cell shape. According to a recent theoretical model of lamellipodia-based dynamics, a similar “two-arc shape” appeared when the protrusive force was high (Sadhu et al., 2023). The dumbbell-shape may thus be a prevalent shape feature that was heretofore overlooked due to peculiarity of the model cells. Indeed, a similar dumbbell-shape has been reported in fragmented Physarum polycephalum (Rieu et al., 2015).

In the E. coli run-and-tumble, the underlying biochemical network has been proposed to be optimally designed to extract binary information in a noisy environment (Nakamura and Kobayashi, 2021). Some bacterial species make use of multiple run modes that differ in how they are modulated in the presence of chemoattractants (Alirezaeizanjani et al., 2020) suggesting diversity and depth at which random walk strategies are likely employed in prokaryotes. In Dictyostelium amoebae, the run length increases under starvation (Haastert and Bosgraaf, 2009) which may be related to their foraging strategy. In immune cells, high correlation between cell speed and persistence is thought to underlie their search efficiency in vivo (Shaebani et al., 2020; Shaebani et al., 2022). Cancer cells show persistent random walk in the metastatic state while weakly persistent in non-metastatic state (Huda et al., 2018). Although chemoattractants for N. gruberi are so far unknown, in Naegleria fowleri, formylated peptides are known to act as chemokine (Marciano-Cabral and Cline, 1987) meaning that it enhances cell polarity and movement in the absence of gradient. Cell-cell variability in such response may explain how a minority of N. gruberi cells under our experimental condition showed persistent monopolarity. Naegleria fowleri are one of several known “brain-eating” amoebae that cause fatal central nervous system infection called amebic meningoencephalitis. Their pathogenicity is thought to be related to their capacity to enter brain by penetrating nasopharygeal mucosa and migrate along olfactory nerves (Thong and Ferrante, 1986). In future works, it should be informative to study how the properties quantified in this work are modulated by chemotactic and chemokinetic factors and how they are related to exploratory and invasive strategies.

4 Materials and methods

4.1 Cell culture

Naegleria gruberi strain NEG-M was obtained from American Type Culture Collection (ATCC 30224). For routine cell propagation, small bits of frozen stock were scraped off using a sharp needle onto a fresh lawn of Klebsilella aerogenes on a NM agar plate (Peptone, Dextrose, K₂HPO₄, KH₂PO4, 2% bactoagar) (Fulton, 1970). The two-member culture plate was incubated at 30°C for a few days until cleared plaques appeared. To start axenic culture, growing cells were picked from the edge of a plaque and suspended in Milli-Q water. 10 μL of the cell suspension was added to 25 mL modified HL5 media (Fulton, 1970) supplemented with 40 ng/mL vitamin B12 and 80 ng/mL folic acid, 10% fetal bovine serum (FBS, Sigma 172,012) and 1% Penicillin-Streptomycin (Gibco) in a 75 cm² canted-neck plastic flask (Corning 431464U). Cells were allowed to attach to the bottom of the flask and incubated at 30°C for 3 days before harvesting for imaging.

4.2 Time-lapse imaging

Axenic growing cells were dislodged from the flask bottom by gentle agitation. Cells were pelleted by centrifugation at 7 × 10² G for 3 min and resuspended in fresh HL5 media. The medium contains 5 mM KH₂PO₄/Na₂HPO buffer and thus provides required electrolytes (King et al., 1979) for optimal migration. The cell density was adjusted to 3.3 × 10² cells/mL for the observations. 3 mL of the cell suspension was plated on a 35 mm glass bottom dish (No. 0 20 mm hole diameter, MatTek). The plate was set to the stage of an inverted microscope (IX81, Olympus) equipped with either a thermal plate or a closed stage-top incubator set to 30°C and left still for 30 min before starting time-lapse image acquisition. All image acquisition was performed at 30°C.

Phase contrast images were obtained by ×40 (LUCPLFLN) objective lens and a sCMOS camera (Prime 95B, Photometrics). To track target cells at multiple non-overlapping fields of view, Micromanager software with a custom written plugin was employed. Timelapse images were obtained from 2 or 3 positions at an interval of 1 s for up to 1 h. Each position was chosen so that initially only a single cell at the center existed in the entire field of view. In between each image acquisition, the cell centroid was calculated from a mask obtained by applying the “Make Binary” function in ImageJ to the most recent image. The automated stage was then recentered to cancel out the centroid displacement.

4.3 Analysis

4.3.1 Characteristics of cellular trajectories

Binary masks from timelapse images were prepared using LABKIT (Arzt et al., 2022). Trajectories of cell centroid were extracted from the mask images using the ImageJ plugin TrackMate (Ershov et al., 2022). The generalized Langevin equation (Eqs 4–6a, b) was numerically solved using the Euler-Maruyama method at 2-milisecond interval with the TorchSDE library (Li et al., 2020; Kidger et al., 2021). Simulated data were sampled at 1 s interval. Velocity $\widehat{\overrightarrow{v}}$ and acceleration $\widehat{\overrightarrow{a}}$ were calculated from the difference in the sampled positions $\widehat{\overrightarrow{r}}$ at time $t_n=n\delta t$ with an interval $\delta t$:

$\widehat{\overrightarrow{v}}\left(t_n\right)=\frac{\left(\widehat{\overrightarrow{r}}\left(t_{n+1}\right)-\widehat{\overrightarrow{r}}\left(t_n\right)\right)}{\delta t}(\mathrm{6a})$

$\widehat{\overrightarrow{a}}\left(t_n\right)=\frac{\left(\widehat{\overrightarrow{v}}\left(t_{n+1}\right)-\widehat{\overrightarrow{v}}\left(t_n\right)\right)}{\delta t}.(\mathrm{6b})$

MSD $msd\left(m\delta t\right)$, probability distribution of speed $p\left(v\right)$, velocity autocorrelation $vac\left(m\mathrm{\delta }t\right)$, mean and standard deviation of acceleration conditional on speed $\left({\widehat{\overline{a}}}_{\parallel }\left(v^{\prime }\right),{\widehat{\overline{a}}}_{\perp }\left(v^{\prime }\right),{\sigma }_{{\widehat{a}}_{\parallel }}\left(v^{\prime }\right),{\sigma }_{{\widehat{a}}_{\perp }}\left(v^{\prime }\right)\right)$, and conditional-averaged strength of turning $〈\cos \theta 〉\left(v^{\prime }\right)$ were calculated from the trajectories for both the experiment and simulation data according to following equations:

$msd\left(m\mathrm{\delta }t\right)={〈{\left(\widehat{\overrightarrow{r}}\left(t_{n+m}\right)-\widehat{\overrightarrow{r}}\left(t_n\right)\right)}^2〉}_n(\mathrm{7a})$

$\kappa \left(v^{\prime }\right)=\left\{\left(i,n\right)|v^{\prime }\le |{\widehat{\overrightarrow{v}}}_i\left(t_n\right)|<v^{\prime }+\delta v\right\}(\mathrm{7b})$

$p\left(v^{\prime }\right)=\frac{\#\kappa \left(v^{\prime }\right)}{\#\left\{\left(i,n\right)\right\}\delta v}(\mathrm{7c})$

$vac\left(m\delta t\right)={〈{\widehat{\overrightarrow{v}}}_i\left(t_n\right)\cdot {\widehat{\overrightarrow{v}}}_i\left(t_{n+m}\right)〉}_{\left(i,n\right)}(\mathrm{7d})$

${\widehat{\overline{a}}}_{\parallel }\left(v^{\prime }\right)={〈{\widehat{a}}_{i,\parallel }\left(t_n\right)〉}_{\left(i,n\right)\in \kappa \left(v^{\prime }\right)},\hspace{1em}{\widehat{\overline{a}}}_{\perp }\left(v^{\prime }\right)={〈{\widehat{a}}_{i,\perp }\left(t_n\right)〉}_{\left(i,n\right)\in \kappa \left(v^{\prime }\right)}(7e)$

${\sigma }_{{\widehat{a}}_{\parallel }}\left(v^{\prime }\right)={std\left({\widehat{a}}_{i,\parallel }\left(t_n\right)\right)}_{\left(i,n\right)\in \kappa \left(v^{\prime }\right)},\hspace{1em}{\sigma }_{{\widehat{a}}_{\perp }}\left(v^{\prime }\right)={std\left({\widehat{a}}_{i,\perp }\left(t_n\right)\right)}_{\left(i,n\right)\in \kappa \left(v^{\prime }\right)}(\mathrm{7f})$

${\widehat{a}}_{i,\parallel }\left(t_n\right)={\widehat{\overrightarrow{a}}}_i\left(t_n\right)\cdot \frac{{\widehat{\overrightarrow{v}}}_i\left(t_n\right)}{\left|{\widehat{\overrightarrow{v}}}_i\left(t_n\right)\right|},\hspace{1em}{\widehat{a}}_{i,\perp }\left(t_n\right)={\widehat{\overrightarrow{a}}}_i\left(t_n\right)\times \frac{{\widehat{\overrightarrow{v}}}_i\left(t_n\right)}{\left|{\widehat{\overrightarrow{v}}}_i\left(t_n\right)\right|}(\mathrm{7g})$

$〈\cos \theta 〉\left(v^{\prime }\right)={\langle \frac{{\widehat{\overrightarrow{v}}}_i\left(t_n\right)}{\left|{\widehat{\overrightarrow{v}}}_i\left(t_n\right)\right|}\cdot \frac{{\widehat{\overrightarrow{v}}}_i\left(t_{n+1}\right)}{\left|{\widehat{\overrightarrow{v}}}_i\left(t_{n+1}\right)\right|}\rangle }_{\left(i,n\right)\in \kappa \left(v^{\prime }\right)}(\mathrm{7h})$

where ${〈〉}_X$ is the average over $X$, subscript $i$ indicates the $i$-th trajectory, $\#$ is the number of items in the following set $\left\{\right\}$, $std$ denotes the unbiased standard deviation, and $\delta v=$ 0.1 μm/s is the bin width. Additionally, we checked the detail of the time evolution of velocity by calculating the autocorrelation of the magnitude and the angle:

$\left|\widehat{\overrightarrow{v}}\right|autocorrelation\left(m\delta t\right)=$

${〈\left(|{\widehat{\overrightarrow{v}}}_i|\left(t_n\right)-{〈|{\widehat{\overrightarrow{v}}}_{i^{\prime }}|\left(t_{n^{\prime }}\right)〉}_{\left(i^{\prime },n^{\prime }\right)}\right)\left(|{\widehat{\overrightarrow{v}}}_i|\left(t_{n+m}\right)-{〈|{\widehat{\overrightarrow{v}}}_{i^{\prime }}|\left(t_{n^{\prime }}\right)〉}_{\left(i^{\prime },n^{\prime }\right)}\right)〉}_{\left(i,n\right)}(\mathrm{8a})$

$\mathit{arg}\left(\widehat{\overrightarrow{v}}\right)autocorrelation\left(m\delta t\right)={\langle \frac{{\widehat{\overrightarrow{v}}}_i\left(t_n\right)}{\left|{\widehat{\overrightarrow{v}}}_i\left(t_n\right)\right|}\cdot \frac{{\widehat{\overrightarrow{v}}}_i\left(t_{n+m}\right)}{\left|{\widehat{\overrightarrow{v}}}_i\left(t_{n+m}\right)\right|}\rangle }_{\left(i,n\right)}.(\mathrm{8b})$

4.3.2 Velocity distribution

We fit a Gaussian distribution to both $v_x$ and $v_y$ to determine the standard deviation ${\sigma }_G$. For $\widehat{\overrightarrow{v}}$ that follows 2-dimensional Gaussian distribution, the distribution of $\left|\widehat{\overrightarrow{v}}\right|$ is readily derived from the chi-square distribution with 2 degrees of freedom where the square of $\left|\widehat{\overrightarrow{v}}\right|/{\sigma }_G$ follows:

$p\left(\frac{{\left|\widehat{\overrightarrow{v}}\right|}^2}{{\sigma }_G^2}\in \left[v^{sq},v^{sq}+d{v}^{sq}\right]\right)=\frac{1}{2}\mathit{exp}\left(-\frac{v^{sq}}{2}\right)d{v}^{sq}\therefore p\left(\left|\widehat{\overrightarrow{v}}\right|\in \left[v^{\prime },v^{\prime }+d{v}^{\prime }\right]\right)=\frac{1}{2}\mathit{exp}\left(-\frac{{v^{\prime }}^2}{2{\sigma }_G^2}\right){\left.\frac{d\left({\left|\widehat{\overrightarrow{v}}\right|}^2/{\sigma }_G^2\right)}{d\left|\widehat{\overrightarrow{v}}\right|}\right|}_{\left|\widehat{\overrightarrow{v}}\right|=v^{\prime }}d{v}^{\prime }$

$=\frac{v^{\prime }}{{\sigma }_G^2}\mathit{exp}\left(-\frac{{v^{\prime }}^2}{2{\sigma }_G^2}\right)d{v}^{\prime }.(9)$

To note, the peak of the above distribution is located at $\left|\widehat{\overrightarrow{v}}\right|={\sigma }_G$.

4.3.3 Fitting VAC

To fit the experimental data with the generalized Langevin equation (Eqs 5a, b), we employed the analytical solution for the velocity autocorrelation ${\mathrm{vac}}^{ss}$. For the observed velocity $\widehat{\overrightarrow{v}}\left(t_n\right)=\delta t^{-1}{\int }_{t_n}^{t_{n+1}}\overrightarrow{v}\left(t\right)dt$, the autocorrelation is:

$va{c}^{ss}\left(m\delta t\right)\equiv {〈\widehat{\overrightarrow{v}}\left(t_{n+m}\right)\cdot \widehat{\overrightarrow{v}}\left(t_n\right)〉}_n^{\text{ss}}={\varphi }_{+}{e}^{-\left(m-1\right){\lambda }_{+}\delta t}{\left(\frac{1-e^{-{\lambda }_{+}\delta t}}{{\lambda }_{+}\delta t}\right)}^2+{\varphi }_{-}{e}^{-\left(m-1\right){\lambda }_{-}\delta t}{\left(\frac{1-e^{-{\lambda }_{-}\delta t}}{{\lambda }_{-}\delta t}\right)}^2(\mathrm{10a})$

${\begin{array}{c}\lambda \end{array}}_{\pm }=\frac{\left(\beta +\gamma \right)\pm \sqrt{{\left(\beta -\gamma \right)}^2+4{\alpha }^2}}{2}(\mathrm{10b})$

${\varphi }_{\pm }={\sigma }^2\frac{1\pm \left(\beta -\gamma \right)}{{\left(\beta -\gamma \right)}^2+4{\alpha }^2}\left(\frac{1\pm \left(\beta -\gamma \right)}{4{\lambda }_{\pm }}+\frac{1\mp \left(\beta -\gamma \right)}{\left(\beta +\gamma \right)}\right).(\mathrm{10c})$

Optimal values of $\alpha ,\beta ,\gamma ,\sigma$ were obtained by minimizing the mean square error between $va{c}^{ss}\left(m\delta t\right)$ and $va{c}^{\mathit{exp}}\left(m\delta t\right)$.

4.3.4 Positional uncertainty

Parameters in Table 1 were obtained by fitting VAC at $\tau \ge 2\mathit{sec}.$ As for the simulation only with generalized Langevin equations (Eqs 5a, b), VAC matched poorly for the shortest time interval of our data ($\tau =$ 0 and 1 s) due to measurement uncertainty arising from finite time step and spatial resolution of the observation. Because acceleration was also defined as the velocity difference in this time interval, the magnitude of acceleration in the simulations was off by one order of magnitude from the real cell data. We emulated these effects in the simulations by including white Gaussian noise with the observed standard deviation ${\sigma }_X$ (see Methods, Table 1). The value of VAC changed only at the shortest time window of $\tau =$ 0 and 1 s by this correction.

To represent positional uncertainty, we incorporated additive noise in the model so that

$\widehat{\overrightarrow{r}}\left(t_n\right)=\overrightarrow{r}\left(t_n\right)+{\sigma }_X{\overrightarrow{\xi }}_n^{\left(X\right)}(11)$

where ${\overrightarrow{\xi }}_n^{\left(X\right)}$ is white gaussian noise which satisfies $〈{\overrightarrow{\xi }}_n^{\left(X\right)}〉=\overrightarrow{0,}〈{\overrightarrow{\xi }}_n^{\left(X\right)}{\overrightarrow{\xi }}_m^{{\left(X\right)}^T}〉=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right){\delta }_{nm}$ where ${\delta }_{nm}$ is the Kronecker delta, and thus independent of all the other variables. ${\sigma }_X$ is the strength of the positional noise, $\overrightarrow{r}\left(t_n\right)$ is the position sampled at time $t_n$, calculated by integrating $\overrightarrow{v}\left(t\right)$ in time according to Eqs 5a, b. The observed velocity ${\widehat{\overrightarrow{v}}}^{\prime }\left(t\right)$ used in the analysis is defined as follows:

${\widehat{\overrightarrow{v}}}^{\prime \left(t_n\right)}=\frac{\widehat{\overrightarrow{r}}\left(t_{n+1}\right)-\widehat{\overrightarrow{r}}\left(t_n\right)}{\delta t}.(12)$

Due to the positional noise, the analytical solution of the velocity autocorrelation at steady state becomes

$va{c}^{ssX}\left(m\delta t\right)\equiv {〈{\widehat{\overrightarrow{v}}}^{\prime }\left(t_{n+m}\right)\cdot {\widehat{\overrightarrow{v}}}^{\prime }\left(t_n\right)〉}_n^{\text{ss}}=\left\{\begin{array}{cc}va{c}^{ss}\left(0\right)+\frac{2{\sigma }_X^2}{\delta t^2}& \left(m=0\right)\\ va{c}^{ss}\left(\mathrm{\delta }t\right)-\frac{{\sigma }_X^2}{\delta t^2}& \left(m=1\right)\\ va{c}^{ss}\left(m\delta t\right)& \left(m\ge 2\right)\end{array}\right..(13)$

The optimal values of $\alpha ,\beta ,\gamma ,\sigma ,{\sigma }_X$ were obtained by minimizing mean square error between $va{c}^{ssX}\left(m\delta t\right)$ and $va{c}^{\mathit{exp}}\left(m\delta t\right)$.

4.3.5 Cell boundary analysis

A MATLAB code for the active contour method (Driscoll et al., 2012)—BoundaryTrack (Nakajima et al., 2016; Fujimori et al., 2019) was used to plot kymographs of the curvature and protrusion velocity of the cell binary mask contour. In brief, the kymographs show time-evolution of curvature or normal vector-projected velocity on the contour. The angle of normal vector was also obtained using this code.

4.3.5.1 Boundary point tracking by BoundaryTrack

Initially, BoundaryTrack detects the sequence of boundary pixels of the mask starting clockwise from the upper-left most point (Supplementary Figure S4A left). At each frame, the boundary was divided into equally spaced 500 points, where the upper-left most point was assigned index 1 (Supplementary Figure S4A center). The boundary points in two consecutive frames were linked so that the mean square of the distance between the linked points was minimized (Supplementary Figure S4A right). As for the latter frame, the index of the point linked with the first point in the previous frame was reset to 1. From the assigned boundary points, the curvature and the velocity were calculated. In particular, the velocity was obtained by calculating the displacement of the points assigned with the same index over time.

4.3.5.2 Comparing the protrusion velocity and the cell centroid velocity

To detect the forward region of the cell, the $i$-th boundary point at time $t$ in the velocity kymograph ${\left\{u_i\left(t\right)\right\}}_{i=1,\dots ,500}$ were smoothed by fitting the velocity values at boundary points in each time with the following joint function:

$\begin{array}{c}{u}_i\left(t\right)=\left\{\begin{array}{cc}{A}_1\left(t\right)\cos \left(\frac{\pi }{L\left\{I\left(t\right)\right\}}\mathit{min}\left(\left|i-C\left\{I\left(t\right)\right\}\right|,\left|i-C\left\{I\left(t\right)\right\}\pm 500\right|\right)\right)& \left(i\in I\left(t\right)\right)\\ -A_2\left(t\right)\cos \left(\frac{\pi }{500-L\left\{I\left(t\right)\right\}}\mathit{min}\left(\left|i-\overline{C}\left\{I\left(t\right)\right\}\right|,\left|i-\overline{C}\left\{I\left(t\right)\right\}\pm 500\right|\right)\right)& \left(i\notin I\left(t\right)\right)\end{array}\right.\hfill \\ I\left(t\right)=\left\{\begin{array}{cc}\left.\left.\left\{i\right.\right|i_1\le i\le i_2\right\}& \left(i_1\le i_2\right)\\ \left\{i\left|\left.i\le i_2\le i_1\right\}\right.\cup \left.\left\{i\right.\right|i_2\le i_1\le i\right\}& \left(i_2<i_1\right)\end{array}\right.\hfill \end{array}$

where $A_1,A_2\ge 0$, $I\left(t\right)$ is a continuous front region bounded by two ends $i_1\left(t\right)$ and $i_2\left(t\right)$. The center $\left\{I\left(t\right)\right\}$, length $L\left\{I\left(t\right)\right\}$, the center of rear region $\overline{C}\left\{I\left(t\right)\right\}$ were defined in the coordinate with the periodic boundary condition, as follows:

$\begin{array}{c}C\left\{I\left(t\right)\right\}=\left\{\begin{array}{cc}\left(i_1+i_2\right)/2& \left(i_1\le i_2\right)\\ \left(i_1+i_2+500\right)/2& \left(i_2<i_1\wedge i_1+i_2\le 500\right)\\ \left(i_1+i_2-500\right)/2& \left(i_2<i_1\wedge i_1+i_2>500\right)\end{array}\right.\overline{C}\left\{I\left(t\right)\right\}=\left\{\begin{array}{cc}\left(i_1+i_2\right)/2& \left(i_2<i_1\right)\\ \left(i_1+i_2+500\right)/2& \left(i_1\le i_2\wedge i_1+i_2\le 500\right)\\ \left(i_1+i_2-500\right)/2& \left(i_1\le i_2\wedge i_1+i_2>500\right)\end{array}\right.\hfill \\ L\left\{I\left(t\right)\right\}=\left\{\begin{array}{cc}{i}_2-i_1& \left(i_1\le i_2\right)\\ 500-\left(i_1-i_2\right)& \left(i_2<i_1\right)\end{array}\right..\hfill \end{array}$

To investigate the relation of front or rear region with the direction of cell centroid velocity, we calculated the angle difference between the normal vector at the center of front or rear region and the centroid velocity.

4.3.5.3 Curvature wave tracking and the leading edge detection

To track the curvature waves, we first detected protrusive regions as follows. Depending on the curvature ${\left\{c_i\left(t\right)\right\}}_{i=1,\dots ,500}$, position #i in the curvature kymograph were classified as either “protrusive” $\left(c_i\left(t\right)>c^{\left(2\right)}\right)$, “flat” $\left(c^{\left(1\right)}<c_i\left(t\right)\le c^{\left(2\right)}\right)$ or “caved” $\left(c_i\left(t\right)\le c^{\left(1\right)}\right)$ where the thresholds $c^{\left(1\right)},c^{\left(2\right)}$ were obtained by the Otsu’s method. At each time point $t$, continuous protrusive regions (j = 1, 2, 3 … ) were defined as set $I_{j,t}^c$ of neighboring protrusive boundary points i between two ends $\left(i_{j,t}^L,i_{j,t}^R\right)\in I^{LR}\left(t\right)$:

$\begin{array}{c}{I}^{LR}\left(t\right)=\left.\left\{\left(i_{j,t}^L,i_{j,t}^R\right)\in {\mathbb{Z}}^2\right.\right|\left(c_{i_{j,t}^L-1}\left(t\right)\le c^{\left(2\right)}\right)\wedge \left(c_{i_{j,t}^R+1}\left(t\right)\le c^{\left(2\right)}\right)\wedge \left(\forall is.t.\left(i_{j,t}^L\le i\le i_{j,t}^R\vee i\le i_{j,t}^R\le i_{j,t}^L\vee i_{j,t}^R\le i_{j,t}^L\le i\right),c_i\left(t\right)>c^{\left(2\right)}\right)\}\hfill \\ I_{j,t}^c=\left\{\begin{array}{cc}\left\{\left.i\right|i_{j,t}^L\le i\le i_{j,t}^R\right\}& \left(i_{j,t}^L\le i_{j,t}^R\right)\\ \left\{i\left|i\le i_{j,t}^R\le i_{j,t}^L\right.\right\}\cup \left\{\left.i\right|i_{j,t}^R\le i_{j,t}^L\le i\right\}& \left(i_{j,t}^R<i_{j,t}^L\right)\end{array}\right.C\left\{I_{j,t}^c\right\}=\left\{\begin{array}{cc}\left(i_{j,t}^L+i_{j,t}^R\right)/2& \left(i_{j,t}^L\le i_{j,t}^R\right)\\ \left(i_{j,t}^L+i_{j,t}^R+500\right)/2& \left(i_{j,t}^R<i_{j,t}^L\wedge i_{j,t}^L+i_{j,t}^R\le 500\right)\\ \left(i_{j,t}^L+i_{j,t}^R-500\right)/2& \left(i_{j,t}^R<i_{j,t}^L\wedge i_{j,t}^L+i_{j,t}^R>500\right)\end{array}\right.\hfill \end{array}$

where $C\left\{I_{j,t}^c\right\}$ is the center of $j$-th protrusive region.

Next, we traced the curvature waves by linking the $j$-th fragment at frame $t$ and the $j^{\prime }$-th fragment at frame $t+1$ if $I_{j,t}^c$ and $I_{j^{\prime },t+1}^c$ have overlapping points. Thus, the set of linked fragments $J^c$ was defined as follows:

$J^c=\left\{\left(j,j^{\prime },t\right)|\exists i\in I_{j^{\prime },t+1}^c,i\in I_{j,t}^c\right\}.$

From each pair of the linked fragments $\left(j,j^{\prime },t\right)\in J^c$, we obtained the angular velocity ${\omega }_{j,j^{\prime }}^c\left(t\right)$ of a protruding region as follows:

${\omega }_{j,j^{\prime }}^c\left(t\right)=\left({\phi }_{C\left\{I_{j^{\prime },t+1}^c\right\}}\left(t+1\right)-{\phi }_{C\left\{I_{j,t}^c\right\}}\left(t\right)\right)/\Delta t$

where ${\phi }_i\left(t\right)$ is the angle of the normal vector at point $i$ at time $t$. The representative angular velocity ${\omega }^c$ were obtained by fitting the histogram of $\left|{\omega }_{j,j^{\prime }}^c\left(t\right)\right|$ to an exponential distribution for all the linked fragments.

To investigate the relation between the curvature wave and the centroid velocity angle, we selected a single dominant wave $j^d\left(t\right)$ whose angle of normal vector ${\phi }_{C\left\{I_{j^d\left(t\right),t}^c\right\}}\left(t\right)$ was closest to that of the centroid velocity at time $t$. The lifetime $\left\{{\tau }_k^d\right\}$ of the leading edge was measured by calculating the time window during which the leading edge was assigned to a particular curvature wave. To this end, we computed the time interval between the time points $t_k^d\in T^d$ at which $j^d\left(t_k^d\right)$ become un-linked to the dominant wave at the next time frame $j^d\left(t_k^d+1\right)$:

$T^d=\left\{t_k^d|\left(j^d\left(t_k^d\right),j^d\left(t_k^d+1\right),t_k^d\right)\notin J^c\right\}{\tau }_k^d=\left(t_{k+1}^d-t_k^d-1\right)\Delta t$

where the index $k$ is given so that $t_k^d$ is listed in the ascending order, i.e., $t_k^d<t_{k+1}^d$ for all integer $k$. We fit a histogram of ${\tau }_k^d$ for all the linked fragments with exponential distribution to obtain the typical duration time of driving wave ${\tau }^d$.

4.3.5.4 Estimating the time scale of centroid velocity autocorrelation

The angular velocity ${\omega }^c$ and the duration time ${\tau }^d$ obtained above were used to estimate the autocorrelation of the angle of cell centroid velocity $\psi \left(t\right)$. The time evolution of $\psi \left(t\right)$ was modeled as 1D persistent random walk with time scale ${\tau }^d$ and step size ${\omega }^c{\tau }^d$. Then the probability distribution of the angle difference $\Delta \psi \in \left(-\infty ,\infty \right)$ is:

$p\left(\psi \left(t\right)-\psi \left(0\right)=\Delta \psi \right)=\frac{1}{\sqrt{\pi {\left({\omega }^c\right)}^2{\tau }^{d}t}}\mathit{exp}\left(-\frac{\Delta {\psi }^2}{{\left({\omega }^c\right)}^2{\tau }^{d}t}\right).$

Therefore, the autocorrelation $\text{AC}\left(\mathrm{t}\right)$ is:

$\text{AC}\left(\mathrm{t}\right)={\int }_{-\infty }^{\infty }p\left(\psi \left(t\right)-\psi \left(0\right)=\Delta \psi \right)\cos \Delta \psi d\Delta \psi =Re\left({\int }_{-\infty }^{\infty }\frac{1}{\sqrt{\pi {\left({\omega }^c\right)}^2{\tau }^{d}t}}\mathit{exp}\left(-\frac{\Delta {\psi }^2}{{\left({\omega }^c\right)}^2{\tau }^{d}t}+i\Delta \psi \right)d\Delta \psi \right)=\mathit{exp}\left(-\frac{{\left({\omega }^c\right)}^2{\tau }^d}{4}t\right).$

Thus, the estimated decay time of the autocorrelation is $T^{est}=4/{\left({\omega }^c\right)}^2{\tau }^d$.

4.3.6 Cell morphology analysis

4.3.6.1 Fourier-based shape analysis

To quantify cell shape, we calculated the elliptic Fourier descriptor (Kuhl and Giardina, 1982). First, we extracted the outline of cell binary mask with a homemade code according to (Nakajima et al., 2016; Fujimori et al., 2019). The periphery of a cell mask $\mathrm{\Gamma }$ was defined as a folded line parametrized with length $0\le \mathcal{l}<L$ connecting the pixels ${\overrightarrow{q}}_i$ on the edge, where each pixel $i$ has pixel $i-1$ and $i+1$ in its 4 nearest neighbor pixels:

$\mathrm{\Gamma }=\left\{\overrightarrow{q}\left(\mathcal{l}\right)|0\le \mathcal{l}<L,\overrightarrow{q}\left(\mathcal{l}\right)=\left\{\begin{array}{c}\begin{array}{cc}{\overrightarrow{q}}_i+\left(\mathcal{l}-\mathrm{i}\right)\left({\overrightarrow{q}}_i-{\overrightarrow{q}}_{i+1}\right)& \left(i\le \mathcal{l}<i+1\right)\end{array}\\ \begin{array}{cc}{\overrightarrow{q}}_{L-1}+\left(\mathcal{l}-\left(L-1\right)\right)\left({\overrightarrow{q}}_{L-1}-{\overrightarrow{q}}_0\right)& \left(L-1\le \mathcal{l}<L\right)\end{array}\end{array}\right\}\right..(14)$

Next, the polygonal outline was converted to 160 equally spaced points ${\left\{{\widehat{\overrightarrow{q}}}_i\right\}}_{i=0}^{159}$ on a relative position on $\mathrm{\Gamma }$:

${\overrightarrow{q}}_i^{\prime }=\sqrt{\frac{3000}{A}}\overrightarrow{q}\left(\frac{i}{160}L\right)(\mathrm{15a})$

${\widehat{\overrightarrow{q}}}_i={\overrightarrow{q}}_i^{\prime }-\frac{1}{160}{\displaystyle \sum _{j=0}^{159}}{\overrightarrow{q}}_j^{\prime }\left(\frac{j}{160}L\right)(\mathrm{15b})$

which is rescaled according to the total number of pixels $A$ in the mask, and the coordinate was set so that the origin is at the cell centroid.

The elliptic Fourier descriptor was calculated by taking the discrete Fourier transformation of ${\widehat{\overrightarrow{r}}}_i$ with wave number $k$:

${\stackrel{\sim }{\overrightarrow{q}}}_k={\displaystyle \sum _{i=0}^{159}}\mathcal{R}\left(-2\pi k\frac{i}{160}\right){\widehat{\overrightarrow{q}}}_i,k=\mathrm{0,1},\dots ,159(16)$

where $\mathcal{R}\left(\cdot \right)$ is a rotational matrix. Its power spectrum $S_k={\left|{\stackrel{\sim }{\overrightarrow{q}}}_k\right|}^2$ was calculated. C_n and C_-n are complex number equivalents of ${\stackrel{\sim }{\overrightarrow{q}}}_{n-1}$ and ${\stackrel{\sim }{\overrightarrow{q}}}_{161-n}$ .

4.3.6.2 Fourier descriptor PCA

We calculated principal component vectors from the representative dataset containing 63 snapshots. From the power spectrum vector $\overrightarrow{S}\equiv \left(S_0,S_1,\dots S_{159}\right)$ for each mask in the representative dataset, averaged power spectrum vector $\overline{\overrightarrow{S}}$ and the covariance matrix $\eta ={\left({\eta }_{kl}\right)}_{k,l=\mathrm{0,1},\dots ,159}$

$\overline{\overrightarrow{S}}=〈\overrightarrow{S}〉(17)$

${\eta }_{kl}=\left\{\begin{array}{c}\begin{array}{cc}var\left(S_k\right)& k=l\end{array}\\ \begin{array}{cc}cov\left(S_k,S_l\right)& k\ne l\end{array}\end{array}\right.(18)$

were calculated, where $var$ and $cov$ denotes the variance and covariance. The $m$-th eigenvalue and eigenvector $\left({\lambda }_m,{\overrightarrow{e}}_m\right)$ of matrix $\eta$ were defined so that the conditions ${\lambda }_1\ge {\lambda }_2\ge \dots \ge {\lambda }_{160}$ and ${\overrightarrow{e}}_m\cdot {\overrightarrow{e}}_{m^{\prime }}={\delta }_{m{m}^{\prime }}$ are met. To note, thus obtained values of $\overline{\overrightarrow{S}}$, ${\lambda }_m$, and ${\overrightarrow{e}}_m$ were used to analyze all the data. Using the eigenvectors, the $m$-th principal component

$PCm=\left(\overrightarrow{S}-\overline{\overrightarrow{S}}\right)\cdot {\overrightarrow{e}}_m(19)$

was calculated for each power spectrum vector of mask.

To characterize cell shape change dynamics, we calculated autocorrelation $A{C}^{PC}$ of PC1 and PC2 values. Using the PC values of cell $i$ at time $t_n$, $A{C}^{PC}$ is:

$\begin{array}{ll}\hfill {\stackrel{\sim }{PC1}}_i\left(t_n\right)& ={PC1}_i\left(t_n\right)-{〈{PC1}_i\left(t_{n^{\prime }}\right)〉}_{n^{\prime }}{\stackrel{\sim }{PC2}}_i\left(t_n\right)={PC2}_i\left(t_n\right)-{〈{PC2}_i\left(t_{n^{\prime }}\right)〉}_{n^{\prime }}A{C}^{PC}\left(m\delta t\right)\hfill \\ \hfill & ={\frac{1}{2}〈{\stackrel{\sim }{PC1}}_i\left(t_n\right){\stackrel{\sim }{PC1}}_i\left(t_{n+m}\right)+{\stackrel{\sim }{PC2}}_i\left(t_n\right){\stackrel{\sim }{PC2}}_i\left(t_{n+m}\right)〉}_{\left(i,n\right)}.\hfill \end{array}$

To restore the shape of cell from a set of principal components $\left(PC1,PC2,\dots ,PC160\right)$, $\overrightarrow{S}$ and ${\widehat{\overrightarrow{q}}}_i$ were sequentially calculated:

$\overrightarrow{S}=\overline{\overrightarrow{S}}+{\displaystyle \sum _{m=1}^{160}}PCm{\overrightarrow{e}}_m(20)$

${\widehat{\overrightarrow{q}}}_i=\frac{1}{160}{\displaystyle \sum _{k=0}^{159}}\mathcal{R}\left(2\pi k\frac{i}{160}\right)\left(\begin{array}{c}\sqrt{S_k}\\ 0\end{array}\right).(21)$

The pixels ${\overrightarrow{q}}_i$ included in the edge were obtained by rounding off ${\widehat{\overrightarrow{q}}}_i$. To show the recovered edge as an image, we made a binary image which has white color only on the pixels ${\overrightarrow{q}}_i$.

4.3.6.3 CNN-based shape analysis

CNN-based PCA and classification were performed based on the morphometry obtained previously (Imoto et al., 2021). In brief, each snapshot image of N. gruberi was input to the pre-trained CNN, and the morphology features were obtained as output. The principal components of these features were calculated using the PCA parameters obtained in (Imoto et al., 2021). The time average of the principal components was taken from all the frames in each time series. According to the morphology features, each snapshot was classified into three morphology classes: Dictyostelium-like, HL60-like, and fish keratocyte-like. Since only two snapshots were classified as keratocyte-like, we conducted the further analysis on Dictyostelium-like, HL60-like classes. The HL60 class ratio was calculated for each timeseries, as the number of snapshots classified as HL60 divided by the total number of snapshots in the timeseries.

4.3.7 Analytical solution of VAC at steady state without positional noise

First, we define VAC as an ensemble-averaged inner product of true velocities at two timepoints:

$vac\left(\Delta t;t\right)=〈\overrightarrow{v}\left(t\right)\cdot \overrightarrow{v}\left(t+\Delta t\right)〉.(22)$

To obtain the dynamics of thus defined VAC, $\overrightarrow{v}\left(t\right)$ can be obtained as itô-integral of generalized Langevin equation with 2-dimensional Brownian motion ${\overrightarrow{B}}_t={\left(B_{x,t},B_{y,t}\right)}^T$:

$d\left(\begin{array}{c}\begin{array}{c}{v}_x\left(t\right)\\ V_x\left(t\right)\end{array}\\ \begin{array}{c}{v}_y\left(t\right)\\ V_y\left(t\right)\end{array}\end{array}\right)=\left(\begin{array}{cc}\begin{array}{cc}-\beta & \alpha \\ \alpha & -\gamma \end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}\\ \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{cc}-\beta & \alpha \\ \alpha & -\gamma \end{array}\end{array}\right)\left(\begin{array}{c}\begin{array}{c}{v}_x\left(t\right)\\ V_x\left(t\right)\end{array}\\ \begin{array}{c}{v}_y\left(t\right)\\ V_y\left(t\right)\end{array}\end{array}\right)dt+\sigma \left(\begin{array}{c}\begin{array}{c}d{B}_{x,t}\\ 0\end{array}\\ \begin{array}{c}d{B}_{y,t}\\ 0\end{array}\end{array}\right)(23)$

$d\left[e^{Ct}\left(\begin{array}{c}{v}_x\left(t\right)\\ V_x\left(t\right)\end{array}\right)\right]=\sigma e^{Ct}\left(\begin{array}{c}d{B}_{x,t}\\ 0\end{array}\right),C=\left(\begin{array}{cc}\beta & -\alpha \\ -\alpha & \gamma \end{array}\right)(\mathrm{24a})$

$d\left[e^{Ct}\left(\begin{array}{c}{v}_y\left(t\right)\\ V_y\left(t\right)\end{array}\right)\right]=\sigma e^{Ct}\left(\begin{array}{c}d{B}_{y,t}\\ 0\end{array}\right)(\mathrm{24b})$

$\therefore \left(\begin{array}{c}{v}_x\left(t\right)\\ V_x\left(t\right)\end{array}\right)=e^{-Ct}\left[\left(\begin{array}{c}{v}_x\left(0\right)\\ V_x\left(0\right)\end{array}\right)+\sigma {\int }_0^t{e}^{C{t}^{\prime }}\left(\begin{array}{c}d{B}_{x,t^{\prime }}\\ 0\end{array}\right)d{t}^{\prime }\right](\mathrm{25a})$

$\left(\begin{array}{c}{v}_y\left(t\right)\\ V_y\left(t\right)\end{array}\right)=e^{-Ct}\left[\left(\begin{array}{c}{v}_y\left(0\right)\\ V_y\left(0\right)\end{array}\right)+\sigma {\int }_0^t{e}^{C{t}^{\prime }}\left(\begin{array}{c}d{B}_{y,t^{\prime }}\\ 0\end{array}\right)d{t}^{\prime }\right].(\mathrm{25b})$

Especially, the velocity can be calculated from the eigenvalues ${\lambda }_{\pm }$ defined above and corresponding eigenvectors ${\overrightarrow{e}}_{\pm }\equiv \left(\begin{array}{c}{e}_{x,\pm }\\ e_{y,\pm }\end{array}\right)$ of $C$ with $e^{Ct}=\left(\begin{array}{cc}{\overrightarrow{e}}_{+}& {\overrightarrow{e}}_{-}\end{array}\right)\left(\begin{array}{cc}{e}^{{\lambda }_{+}t}& 0\\ 0& e^{{\lambda }_{-}t}\end{array}\right)\left(\begin{array}{c}{\overrightarrow{e}}_{+}^T\\ {\overrightarrow{e}}_{-}^T\end{array}\right)$:

$v_x\left(t\right)={\left(\begin{array}{c}{e}_{x,+}\\ e_{x,-}\end{array}\right)}^T\left[\left(\begin{array}{c}{e}^{-{\lambda }_{+}t}{\overrightarrow{e}}_{+}^T\\ e^{-{\lambda }_{-}t}{\overrightarrow{e}}_{-}^T\end{array}\right)\left(\begin{array}{c}{v}_x\left(0\right)\\ V_x\left(0\right)\end{array}\right)+\sigma {\int }_0^t\left(\begin{array}{c}{e}^{{\lambda }_{+}\left(t^{\prime }-t\right)}{e}_{x,+}\\ e^{{\lambda }_{-}\left(t^{\prime }-t\right)}{e}_{x,-}\end{array}\right)d{B}_{x,t^{\prime }}\right](\mathrm{26a})$

$v_y\left(t\right)={\left(\begin{array}{c}{e}_{x,+}\\ e_{x,-}\end{array}\right)}^T\left[\left(\begin{array}{c}{e}^{-{\lambda }_{+}t}{\overrightarrow{e}}_{+}^T\\ e^{-{\lambda }_{-}t}{\overrightarrow{e}}_{-}^T\end{array}\right)\left(\begin{array}{c}{v}_y\left(0\right)\\ V_y\left(0\right)\end{array}\right)+\sigma {\int }_0^t\left(\begin{array}{c}{e}^{{\lambda }_{+}\left(t^{\prime }-t\right)}{e}_{x,+}\\ e^{{\lambda }_{-}\left(t^{\prime }-t\right)}{e}_{x,-}\end{array}\right)d{B}_{y,t^{\prime }}\right](\mathrm{26b})$

${\overrightarrow{e}}_{\pm }=\left(\begin{array}{c}\cos {\theta }_{\pm }\\ \sin {\theta }_{\pm }\end{array}\right),\tan {\theta }_{\pm }=\frac{\beta -\gamma \mp \sqrt{{\left(\beta -\gamma \right)}^2+4{\alpha }^2}}{2\alpha }.(\mathrm{26c})$

Using the representation of $\overrightarrow{v}\left(t\right)$ above and the property of Brownian motion $\int f\left(t\right)d{\overrightarrow{B}}_t=\overrightarrow{0},d{\overrightarrow{B}}_{t}d{\overrightarrow{B}}_{t^{\prime }}^T=\delta \left(t-t^{\prime }\right)\left(\begin{array}{cc}dt& 0\\ 0& dt\end{array}\right)$, VAC is:

$vac\left(\Delta t;t_0\right)=$

$\langle {\displaystyle \sum _{\chi =x,y}}\left\{\left[{\left(\begin{array}{c}{e}_{x,+}\\ e_{x,-}\end{array}\right)}^T\left(\begin{array}{c}{e}^{-{\lambda }_{+}{t}_0}{\overrightarrow{e}}_{+}^T\\ e^{-{\lambda }_{-}{t}_0}{\overrightarrow{e}}_{-}^T\end{array}\right)\left(\begin{array}{c}{v}_{\chi }\left(0\right)\\ V_{\chi }\left(0\right)\end{array}\right)\right]\left[{\left(\begin{array}{c}{e}_{x,+}\\ e_{x,-}\end{array}\right)}^T\left(\begin{array}{c}{e}^{-{\lambda }_{+}\left(t_0+\Delta t\right)}{\overrightarrow{e}}_{+}^T\\ e^{-{\lambda }_{-}\left(t_0+\Delta t\right)}{\overrightarrow{e}}_{-}^T\end{array}\right)\left(\begin{array}{c}{v}_{\chi }\left(0\right)\\ V_{\chi }\left(0\right)\end{array}\right)\right]\right\}$

$+2{\sigma }^2{\int }_0^{t_0}\left[{\left(\begin{array}{c}{e}_{x,+}\\ e_{x,-}\end{array}\right)}^T\left(\begin{array}{c}{e}^{{\lambda }_{+}\left(t^{\prime }-t_0\right)}{e}_{x,+}\\ e^{{\lambda }_{-}\left(t^{\prime }-t_0\right)}{e}_{x,-}\end{array}\right)\right]\left[{\left(\begin{array}{c}{e}_{x,+}\\ e_{x,-}\end{array}\right)}^T\left(\begin{array}{c}{e}^{{\lambda }_{+}\left(t^{\prime }-\left(t_0+\Delta t\right)\right)}{e}_{x,+}\\ e^{{\lambda }_{-}\left(t^{\prime }-\left(t_0+\Delta t\right)\right)}{e}_{x,-}\end{array}\right)\right]d{t}^{\prime }\rangle .(27)$

Since ${\lambda }_{\pm }$ is always positive when $\alpha ,\beta ,\gamma >0$, the first term of VAC disappears with time at the rate of $e^{-2{\lambda }_{-}{t}_0}$. The lower limit of integration also disappears at the same rate. The only time-independent term comes from the upper limit of integration and is the steady state solution of VAC:

$vac\left(\Delta t;t\right){\displaystyle \stackrel{t\to \infty }{\to }}{\varphi }_{+}{e}^{-{\lambda }_{+}\Delta t}+{\varphi }_{-}{e}^{-{\lambda }_{-}\Delta t}(\mathrm{28a})$

${\varphi }_{\pm }=2{\sigma }^2\left(\frac{e_{x,\pm }^4}{2{\lambda }_{\pm }}+\frac{e_{x,+}^2{e}_{x,-}^2}{{\lambda }_{+}+{\lambda }_{-}}\right)(\mathrm{28b})$

where the second line is another representation of ${\varphi }_{\pm }$ defined above. Finally, considering the sampling procedure where the velocity is observed as the difference of discretely sampled positions, the representation of $va{c}^{ss}$ is obtained by time integration of VAC:

$\begin{array}{c}va{c}^{ss}\left(k\delta t\right)=\frac{1}{\delta t^2}{〈\left(\overrightarrow{r}\left(t_{n+k+1}\right)-\overrightarrow{r}\left(t_{n+k}\right)\right)\cdot \left(\overrightarrow{r}\left(t_{n+1}\right)-\overrightarrow{r}\left(t_n\right)\right)〉}_n=\frac{1}{\delta t^2}{〈\left({\int }_{t_{n+k}}^{t_{n+k+1}}\overrightarrow{v}\left(t^{\prime }\right)d{t}^{\prime }\right)\cdot \left({\int }_{t_n}^{t_{n+1}}\overrightarrow{v}\left(t^{″}\right)d{t}^{″}\right)〉}_n\\ =\frac{1}{\delta t^2}{\int }_{k\delta t}^{\left(k+1\right)\delta t}d{t}^{\prime }{\int }_0^{\delta t}d{t}^{″}{〈vac\left(t^{\prime }-t^{″};t_n\right)〉}_n={\varphi }_{+}{e}^{-\left(k-1\right){\lambda }_{+}\delta t}{\left(\frac{1-e^{-{\lambda }_{+}\delta t}}{{\lambda }_{+}\delta t}\right)}^2+{\varphi }_{-}{e}^{-\left(k-1\right){\lambda }_{-}\delta t}{\left(\frac{1-e^{-{\lambda }_{-}\delta t}}{{\lambda }_{-}\delta t}\right)}^2.\end{array}(29)$

In the third raw, we used the relation ${〈vac\left(t^{\prime }-t^{″};t_n\right)〉}_n=\underset{n\to \infty }{\lim }vac\left(t^{\prime }-t^{″};t_n\right)=\underset{t\to \infty }{\lim }vac\left(t^{\prime }-t^{″};t\right)$ because time average should converge to the steady state solution if the VAC itself converges.

Data availability statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author contributions

MU: Data curation, Formal Analysis, Investigation, Methodology, Software, Supervision, Validation, Visualization, Writing–original draft, Funding acquisition. YM: Data curation, Formal Analysis, Investigation, Writing–review and editing, Software. AK: Data curation, Writing–review and editing. DI: Writing–review and editing, Data curation, Investigation, Software, Visualization. SS: Conceptualization, Funding acquisition, Methodology, Project administration, Resources, Supervision, Writing–original draft.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This was work was supported by JSPS KAKENHI Grant Number JP19H05801 to SS, JP22H05673 to MU, JST CREST JPMJCR1923 to SS and partly by HFSP Research Grant RGP0051/2021 to SS.

Acknowledgments

The authors thank the present and the past members of the SS lab for experimental supports and discussion.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2023.1274127/full#supplementary-material

References

Alirezaeizanjani, Z., Großmann, R., Pfeifer, V., Hintsche, M., and Beta, C. (2020). Chemotaxis strategies of bacteria with multiple run modes. Sci. Adv. 6, eaaz6153. doi:10.1126/sciadv.aaz6153

PubMed Abstract | CrossRef Full Text | Google Scholar

Ariel, G., Rabani, A., Benisty, S., Partridge, J. D., Harshey, R. M., and Be’er, A. (2015). Swarming bacteria migrate by Lévy Walk. Nat. Commun. 6, 8396. doi:10.1038/ncomms9396

PubMed Abstract | CrossRef Full Text | Google Scholar

Arzt, M., Deschamps, J., Schmied, C., Pietzsch, T., Schmidt, D., Tomancak, P., et al. (2022). LABKIT: labeling and segmentation toolkit for big image data. Front. Comput. Sci. 4. doi:10.3389/fcomp.2022.777728

CrossRef Full Text | Google Scholar

Bartumeus, F., Catalan, J., Fulco, U. L., Lyra, M. L., and Viswanathan, G. M. (2002). Optimizing the encounter rate in biological interactions: Lévy versus brownian strategies. Phys. Rev. Lett. 88, 097901. doi:10.1103/physrevlett.88.097901

PubMed Abstract | CrossRef Full Text | Google Scholar

Brown, M. W., Silberman, J. D., and Spiegel, F. W. (2012). A contemporary evaluation of the acrasids (Acrasidae, Heterolobosea, Excavata). Eur. J. Protistol. 48, 103–123. doi:10.1016/j.ejop.2011.10.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Chi, Q., Yin, T., Gregersen, H., Deng, X., Fan, Y., Zhao, J., et al. (2014). Rear actomyosin contractility-driven directional cell migration in three-dimensional matrices: a mechano-chemical coupling mechanism. J. Roy. Soc. Interface 11, 20131072. doi:10.1098/rsif.2013.1072

PubMed Abstract | CrossRef Full Text | Google Scholar

Dieterich, P., Klages, R., Preuss, R., and Schwab, A. (2008). Anomalous dynamics of cell migration. Proc. Natl. Acad. Sci. U. S. A. 105, 459–463. doi:10.1073/pnas.0707603105

PubMed Abstract | CrossRef Full Text | Google Scholar

Driscoll, M. K., McCann, C., Kopace, R., Homan, T., Fourkas, J. T., Parent, C., et al. (2012). Cell shape dynamics: from waves to migration. Plos Comput. Biol. 8, e1002392. doi:10.1371/journal.pcbi.1002392

PubMed Abstract | CrossRef Full Text | Google Scholar

Dunn, G. A., and Brown, A. F. (1987). A unified approach to analysing cell motility. J. Cell Sci. 8, 81–102. doi:10.1242/jcs.1987.supplement_8.5

PubMed Abstract | CrossRef Full Text | Google Scholar

Ebata, H., Yamamoto, A., Tsuji, Y., Sasaki, S., Moriyama, K., Kuboki, T., et al. (2018). Persistent random deformation model of cells crawling on a gel surface. Sci. Rep. 8, 5153. doi:10.1038/s41598-018-23540-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Ershov, D., Phan, M.-S., Pylvänäinen, J. W., Rigaud, S. U., Blanc, L. L., Charles-Orszag, A., et al. (2022). TrackMate 7: integrating state-of-the-art segmentation algorithms into tracking pipelines. Nat. Methods 19, 829–832. doi:10.1038/s41592-022-01507-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Fritz-Laylin, L. K., Lord, S. J., and Mullins, R. D. (2017). WASP and SCAR are evolutionarily conserved in actin-filled pseudopod-based motility. J. Cell Biol. 216, 1673–1688. doi:10.1083/jcb.201701074

PubMed Abstract | CrossRef Full Text | Google Scholar

Fritz-Laylin, L. K., Prochnik, S. E., Ginger, M. L., Dacks, J. B., Carpenter, M. L., Field, M. C., et al. (2010). The genome of Naegleria gruberi illuminates early eukaryotic versatility. Cell 140, 631–642. doi:10.1016/j.cell.2010.01.032

PubMed Abstract | CrossRef Full Text | Google Scholar

Fujimori, T., Nakajima, A., Shimada, N., and Sawai, S. (2019). Tissue self-organization based on collective cell migration by contact activation of locomotion and chemotaxis. Proc. Natl. Acad. Sci. U. S. A. 116, 4291–4296. doi:10.1073/pnas.1815063116

PubMed Abstract | CrossRef Full Text | Google Scholar

Fulton, C. (1970). Amebo-flagellates as research partners: the laboratory biology of Naegleria and Tetramitus. Methods Cell Physiology 4, 341–476. doi:10.1016/S0091-679X(08)61759-8

CrossRef Full Text | Google Scholar

Gail, M. H., and Boone, C. W. (1970). The locomotion of mouse fibroblasts in tissue culture. Biophys. J. 10, 980–993. doi:10.1016/s0006-3495(70)86347-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Haastert, P. J. M. V. (2021). Short- and long-term memory of moving amoeboid cells. Plos One 16, e0246345. doi:10.1371/journal.pone.0246345

PubMed Abstract | CrossRef Full Text | Google Scholar

Haastert, P. J. M. V., and Bosgraaf, L. (2009). Food searching strategy of amoeboid cells by starvation induced run length extension. Plos One 4, e6814. doi:10.1371/journal.pone.0006814

PubMed Abstract | CrossRef Full Text | Google Scholar

Harris, T. H., Banigan, E. J., Christian, D. A., Konradt, C., Wojno, E. D. T., Norose, K., et al. (2012). Generalized Lévy walks and the role of chemokines in migration of effector CD8+ T cells. Nature 486, 545–548. doi:10.1038/nature11098

PubMed Abstract | CrossRef Full Text | Google Scholar

Hartman, R., Lau, K., Chou, W., and Coates, T. D. (1994). The fundamental motor of the human neutrophil is not random: evidence for local non-Markov movement in neutrophils. Biophys. J. 67, 2535–2545. doi:10.1016/S0006-3495(94)80743-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, C. -H., Tang, M., Shi, C., Lglesias, P. A., and Devreotes, P. N. (2013). An excitable signal integrator couples to an idling cytoskeletal oscillator to drive cell migration. Nat. Cell Biol. 15, 1307–1316. doi:10.1038/ncb2859

PubMed Abstract | CrossRef Full Text | Google Scholar

Huda, S., Weigelin, B., Wolf, K., Tretiakov, K. V., Polev, K., Wilk, G., et al. (2018). Lévy-like movement patterns of metastatic cancer cells revealed in microfabricated systems and implicated in vivo. Nat. Commun. 9, 4539. doi:10.1038/s41467-018-06563-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Huo, H., He, R., Zhang, R., and Yuan, J. (2021). Swimming Escherichia coli cells explore the environment by Lévy walk. Appl. Environ. Microb. 87, 024299–e2520. doi:10.1128/aem.02429-20

CrossRef Full Text | Google Scholar

Imoto, D., Saito, N., Nakajima, A., Honda, G., Ishida, M., Sugita, T., et al. (2021). Comparative mapping of crawling-cell morphodynamics in deep learning-based feature space. Plos Comput. Biol. 17, e1009237. doi:10.1371/journal.pcbi.1009237

PubMed Abstract | CrossRef Full Text | Google Scholar

Itakura, A., Aslan, J. E., Kusanto, B. T., Phillips, K. G., Porter, J. E., Newton, P. K., et al. (2013). p21-Activated kinase (PAK) regulates cytoskeletal reorganization and directional migration in human neutrophils. Plos One 8, e73063. doi:10.1371/journal.pone.0073063

PubMed Abstract | CrossRef Full Text | Google Scholar

Jerison, E. R., and Quake, S. R. (2020). Heterogeneous T cell motility behaviors emerge from a coupling between speed and turning in vivo. Elife 9, e53933. doi:10.7554/elife.53933

PubMed Abstract | CrossRef Full Text | Google Scholar

Kidger, P., Foster, J., Li, X., Oberhauser, H., and Lyons, T. (2021). Neural SDEs as infinite-dimensional GANs. Available from: https://arxiv.org/abs/2102.03657.

Google Scholar

King, C. A., Davies, A. H., and Preston, T. M. (1981). Lack of substrate specificity on the speed of amoeboid locomotion inNaegleria gruberi. Experientia 37, 709–710. doi:10.1007/bf01967936

CrossRef Full Text | Google Scholar

King, C. A., Westwood, R., Cooper, L., and Preston, T. M. (1979). Speed of locomotion of the soil amoebaNaegleria gruberi in media of different ionic compositions with special reference to interactions with the substratum. Protoplasma 99, 323–334. doi:10.1007/bf01275804

CrossRef Full Text | Google Scholar

Kuhl, F. P., and Giardina, C. R. (1982). Elliptic Fourier features of a closed contour. Comput. Vis. Graph 18, 236–258. doi:10.1016/0146-664x(82)90034-x

CrossRef Full Text | Google Scholar

Kwon, T., Kwon, O.-S., Cha, H.-J., and Sung, B. J. (2019). Stochastic and heterogeneous cancer cell migration: experiment and theory. Sci. Rep. 9, 16297. doi:10.1038/s41598-019-52480-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, L., Cox, E. C., and Flyvbjerg, H. (2011). “Dicty dynamics”: dictyostelium motility as persistent random motion. Phys. Biol. 8, 046006. doi:10.1088/1478-3975/8/4/046006

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, L., Nørrelykke, S. F., and Cox, E. C. (2008). Persistent cell motion in the absence of external signals: a search strategy for eukaryotic cells. Plos One 3, e2093. doi:10.1371/journal.pone.0002093

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X., Wong, T.-K. L., Chen, R. T. Q., and Duvenaud, D. (2020). Scalable gradients for stochastic differential equations. Available from: https://arxiv.org/abs/2001.01328.

Google Scholar

Marciano-Cabral, F., and Cline, M. (1987). Chemotaxis by Naegleria fowleri for bacteria. J. Protozool. 34, 127–131. doi:10.1111/j.1550-7408.1987.tb03147.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Morales, J. C. F., Xue, Q., and Roh-Johnson, M. (2022). An evolutionary and physiological perspective on cell-substrate adhesion machinery for cell migration. Front. Cell Dev. Biol. 10, 943606. doi:10.3389/fcell.2022.943606

PubMed Abstract | CrossRef Full Text | Google Scholar

Nakajima, A., Ishida, M., Fujimori, T., Wakamoto, Y., and Sawai, S. (2016). The microfluidic lighthouse: an omnidirectional gradient generator. Lab. Chip 16, 4382–4394. doi:10.1039/c6lc00898d

PubMed Abstract | CrossRef Full Text | Google Scholar

Nakamura, K., and Kobayashi, T. J. (2021). Connection between the bacterial chemotactic network and optimal filtering. Phys. Rev. Lett. 126, 128102. doi:10.1103/physrevlett.126.128102

PubMed Abstract | CrossRef Full Text | Google Scholar

Ohta, T., Tarama, M., and Sano, M. (2016). Simple model of cell crawling. Physica D. 318–319, 3–11. doi:10.1016/j.physd.2015.10.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Parfrey, L. W., Lahr, D. J. G., Knoll, A. H., and Katz, L. A. (2011). Estimating the timing of early eukaryotic diversification with multigene molecular clocks. Proc. Natl. Acad. Sci. U. S. A. 108, 13624–13629. doi:10.1073/pnas.1110633108

PubMed Abstract | CrossRef Full Text | Google Scholar

Partridge, B., Anton, S. G., Khorshed, R., Adams, G., Pospori, C., Celso, C. L., et al. (2022). Heterogeneous run-and-tumble motion accounts for transient non-Gaussian super-diffusion in haematopoietic multi-potent progenitor cells. Plos One 17, e0272587. doi:10.1371/journal.pone.0272587

PubMed Abstract | CrossRef Full Text | Google Scholar

Passucci, G., Brasch, M. E., Henderson, J. H., Zaburdaev, V., and Manning, M. L. (2019). Identifying the mechanism for superdiffusivity in mouse fibroblast motility. Plos Comput. Biol. 15, e1006732. doi:10.1371/journal.pcbi.1006732

PubMed Abstract | CrossRef Full Text | Google Scholar

Pollard, T. D. (2007). Regulation of actin filament assembly by Arp2/3 complex and formins. Annu. Rev. Bioph. Biom. 36, 451–477. doi:10.1146/annurev.biophys.35.040405.101936

PubMed Abstract | CrossRef Full Text | Google Scholar

Preston, T. M., and King, C. A. (1978). An experimental study of the interaction between the soil amoeba Naegleria gruberi and a glass substrate during amoeboid locomotion. J. Cell Sci. 34, 145–158. doi:10.1242/jcs.34.1.145

PubMed Abstract | CrossRef Full Text | Google Scholar

Prostak, S. M., Robinson, K. A., Titus, M. A., and Fritz-Laylin, L. K. (2021). The actin networks of chytrid fungi reveal evolutionary loss of cytoskeletal complexity in the fungal kingdom. Curr. Biol. 31, 1192–1205.e6. doi:10.1016/j.cub.2021.01.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Reynolds, A. M. (2010). Can spontaneous cell movements be modelled as Lévy walks? Phys. Stat. Mech. Appl. 389, 273–277. doi:10.1016/j.physa.2009.09.027

CrossRef Full Text | Google Scholar

Reynolds, A. M., and Ouellette, N. T. (2016). Swarm dynamics may give rise to Lévy flights. Sci. Rep. 6, 30515. doi:10.1038/srep30515

PubMed Abstract | CrossRef Full Text | Google Scholar

Rieu, J.-P., Delanoë-Ayari, H., Takagi, S., Tanaka, Y., and Nakagaki, T. (2015). Periodic traction in migrating large amoeba of Physarum polycephalum. J. Roy. Soc. Interface 12, 20150099. doi:10.1098/rsif.2015.0099

PubMed Abstract | CrossRef Full Text | Google Scholar

Rodiek, B., and Hauser, M. J. B. (2015). Migratory behaviour of Physarum polycephalum microplasmodia. Eur. Phys. J. Spec. Top. 224, 1199–1214. doi:10.1140/epjst/e2015-02455-2

CrossRef Full Text | Google Scholar

Sadhu, R. K., Iglič, A., and Gov, N. S. (2023). A minimal cell model for lamellipodia-based cellular dynamics and migration. J. Cell Sci. 136, jcs260744. doi:10.1242/jcs.260744

PubMed Abstract | CrossRef Full Text | Google Scholar

Schwarz, G. (1978). Estimating the dimension of a model. Ann. Stat. 6. doi:10.1214/aos/1176344136

CrossRef Full Text | Google Scholar

Sebé-Pedrós, A., Grau-Bové, X., Richards, T. A., and Ruiz-Trillo, I. (2014). Evolution and classification of myosins, a paneukaryotic whole-genome approach. Genome Biol. Evol. 6, 290–305. doi:10.1093/gbe/evu013

PubMed Abstract | CrossRef Full Text | Google Scholar

Selmeczi, D., Li, L., Pedersen, L. I. I., Nrrelykke, S. F., Hagedorn, P. H., Mosler, S., et al. (2008). Cell motility as random motion: a review. Eur. Phys. J. Spec. Top. 157, 1–15. doi:10.1140/epjst/e2008-00626-x

CrossRef Full Text | Google Scholar

Selmeczi, D., Mosler, S., Hagedorn, P. H., Larsen, N. B., and Flyvbjerg, H. (2005). Cell motility as persistent random motion: theories from experiments. Biophys. J. 89, 912–931. doi:10.1529/biophysj.105.061150

PubMed Abstract | CrossRef Full Text | Google Scholar

Shaebani, M. R., Jose, R., Santen, L., Stankevicins, L., and Lautenschläger, F. (2020). Persistence-speed coupling enhances the search efficiency of migrating immune cells. Phys. Rev. Lett. 125, 268102. doi:10.1103/physrevlett.125.268102

PubMed Abstract | CrossRef Full Text | Google Scholar

Shaebani, M. R., Piel, M., and Lautenschläger, F. (2022). Distinct speed and direction memories of migrating dendritic cells diversify their search strategies. Biophys. J. 121, 4099–4108. doi:10.1016/j.bpj.2022.09.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Stokes, C. L., Lauffenburger, D. A., and Williams, S. K. (1991). Migration of individual microvessel endothelial cells: stochastic model and parameter measurement. J. Cell Sci. 99, 419–430. doi:10.1242/jcs.99.2.419

PubMed Abstract | CrossRef Full Text | Google Scholar

Takagi, H., Sato, M. J., Yanagida, T., and Ueda, M. (2008). Functional analysis of spontaneous cell movement under different physiological conditions. Plos One 3, e2648. doi:10.1371/journal.pone.0002648

PubMed Abstract | CrossRef Full Text | Google Scholar

Taktikos, J., Stark, H., and Zaburdaev, V. (2013). How the motility pattern of bacteria affects their dispersal and chemotaxis. Plos One 8, e81936. doi:10.1371/journal.pone.0081936

PubMed Abstract | CrossRef Full Text | Google Scholar

Thong, Y. H., and Ferrante, A. (1986). Migration patterns of pathogenic and nonpathogenic Naegleria spp. Infect. Immun. 51, 177–180. doi:10.1128/IAI.51.1.177-180.1986

PubMed Abstract | CrossRef Full Text | Google Scholar

Tweedy, L., Meier, B., Stephan, J., Heinrich, D., and Endres, R. G. (2013). Distinct cell shapes determine accurate chemotaxis. Sci. Rep. 3, 2606. doi:10.1038/srep02606

PubMed Abstract | CrossRef Full Text | Google Scholar

Velle, K. B., and Fritz-Laylin, L. K. (2020). Conserved actin machinery drives microtubule-independent motility and phagocytosis in Naegleria. J. Cell Biol. 219, e202007158. doi:10.1083/jcb.202007158

PubMed Abstract | CrossRef Full Text | Google Scholar

Velle, K. B., and Fritz-Laylin, L. K. (2019). Diversity and evolution of actin-dependent phenotypes. Curr. Opin. Genet. Dev. 58, 40–48. doi:10.1016/j.gde.2019.07.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Viswanathan, G. M., Buldyrev, S. V., Havlin, S., Luz, M. G. E. D., Raposo, E. P., and Stanley, H. E. (1999). Optimizing the success of random searches. Nature 401, 911–914. doi:10.1038/44831

PubMed Abstract | CrossRef Full Text | Google Scholar

Walsh, C. J. (2007). The role of actin, actomyosin and microtubules in defining cell shape during the differentiation of Naegleria amebae into flagellates. Eur. J. Cell Biol. 86, 85–98. doi:10.1016/j.ejcb.2006.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Welf, E. S., Ahmed, S., Johnson, H. E., Melvin, A. T., and Haugh, J. M. (2012). Migrating fibroblasts reorient directionality by a metastable, PI3K-dependent mechanism. J. Cell Biol. 197, 105–114. doi:10.1083/jcb.201108152

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, X., Wen, X., Veltman, D. M., Keizer-Gunnink, I., Pots, H., Kortholt, A., et al. (2017). GPCR-controlled membrane recruitment of negative regulator C2GAP1 locally inhibits Ras signaling for adaptation and long-range chemotaxis. Proc. Natl. Acad. Sci. U. S. A. 114, E10092–E10101. doi:10.1073/pnas.1703208114

PubMed Abstract | CrossRef Full Text | Google Scholar