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).

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.

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).

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 v → : v → t + τ ∙ v → t t = Φ 1 e − τ T 1 + Φ 2 e − τ T 2 (2) we obtained decay time T ₁ and T ₂ of approximately 6 s and 90 s, respectively, where the weight Φ 1 and Φ 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 v → in 1 s interval) (Supplementary Figures S2J, K). The mean acceleration orthogonal ( a ⊥ ) 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 ∥ ) 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.

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) d v → t d t = − β v → t + α 2 ∫ − ∞ t e − γ t − t ′ v → t ′ d t ′ + σ ξ → t (3)

Here, α is the strength of memory effect, and ξ → t is a normalized Gaussian white noise that satisfies ξ → t = 0 , ξ → t ξ → t ′ T = 1 0 0 1 δ t − t ′ ,   is an ensemble average and δ t is the delta function, σ is the noise strength (Selmeczi et al., 2005). By introducing V → t = α ∫ − ∞ t e − γ t − t ′ v → t ′ d t ′ (4) the equation of motion becomes d v → t d t = − β v → t + α V → t + σ ξ → t (5a) d V → t d t = α v → t − γ V → t (5b)

Based on the values of T 1 , T 2 , Φ 1 , Φ 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 v → ) and the orientation ( v → / v → ) (Supplementary Figure S3). In the experimental data, it is only the autocorrelation of the orientation v → / v → not v → 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 v → / v → and v → decayed with the two exponents (Supplementary Figures S3C, D).

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).

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).

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.

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 PC 1 fourier / d t but not with d PC 2 fourier / d t (Supplementary Figures S7C, D). As the former relation was seen in the negative direction d PC 1 fourier / d t <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 n m ≡ C ˙ − n C m − C − n 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 = C ˙ − 2 C 3 − C − 2 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 = C ˙ − 3 C 4 − C − 3 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.

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.

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.

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.