Systemic cellular migration: The forces driving the directed locomotion movement of cells

Abstract Directional motility is an essential property of cells. Despite its enormous relevance in many fundamental physiological and pathological processes, how cells control their locomotion movements remains an unresolved question. Here, we have addressed the systemic processes driving the directed locomotion of cells. Specifically, we have performed an exhaustive study analyzing the trajectories of 700 individual cells belonging to three different species (Amoeba proteus, Metamoeba leningradensis, and Amoeba borokensis) in four different scenarios: in absence of stimuli, under an electric field (galvanotaxis), in a chemotactic gradient (chemotaxis), and under simultaneous galvanotactic and chemotactic stimuli. All movements were analyzed using advanced quantitative tools. The results show that the trajectories are mainly characterized by coherent integrative responses that operate at the global cellular scale. These systemic migratory movements depend on the cooperative nonlinear interaction of most, if not all, molecular components of cells.


Introduction
Self-locomotion is one of the most important complex behaviors of cells endowed with migratory responses.In the permanent struggle for survival, free cells move efficiently to find food following adequate direction and speed, avoiding predators and adverse conditions.Cell motility is crucial for life in Metazoan organisms and fundamental to establish the appropriate organization of all multicellular organisms, playing a central role in a plethora of essential biological phenomena such as embryogenesis, morphogenesis, organogenesis, neural development, adult tissue remodeling, wound healing, immune responses, angiogenesis, tissue regeneration and repair, cell differentiation, etc. (1).Moreover, the deregulation of cell movements in humans is involved in many pathological processes such as metastatic tumor progression (2)(3)(4)(5), atherosclerosis and other vascular diseases (6), congenital brain pathologies (7), osteoarthritis (8), rheumatoid arthritis (8,9), hearing disorders (10), asthma (11)(12)(13), chronic obstructive pulmonary disease (14), multiple sclerosis (15,16), psoriasis (17,18), Crohn's disease (19,20), and immune-related actinopathies (21).Although cell locomotion was already observed in 1675 by van Leeuwenhoek in his respected pioneer microscopic studies (22), researchers and scholars have not yet come to unveil how cells migrate in the presence of complex cues.
Given its importance, great attention has been focused on the study of the diverse molecular parts involved in directional motility.A relevant number of these experimental studies have unequivocally shown that locomotion movements are complex processes that involve practically all cellular components.So, directed movements are primarily driven by the cytoskeleton (the essential part of the locomotion system) which is a sophisticated dynamic structure formed by three main molecular components: actin microfilaments, microtubules, and intermediate filaments, all of them interacting in complex dynamic networks (23).
In particular, the activity of actin cytoskeleton networks is largely dependent on a wide variety of regulatory molecules such as small GTPases (24), integrins (25), and many posttranslational modifications such as phosphorylation, acetylation, arginylation, oxidation, and others (26).In addition, the dynamic turnover of the actin filament networks is essential to regulate cell migration (27).Cytoskeleton networks are coupled with other complex regulated systems such as membrane surface receptors and signal transduction pathways which also participate in the control of locomotion movements (28).Energy is another essential element in cell motility; when cells move the cytoskeleton transforms chemical energy into mechanical forces (dynein cytoskeletal motor proteins) entailing considerable bioenergetic demands; for such a purpose the mitochondrial activity and the adenylate energy system are important regulators of directionality motion (29).Cell membrane activities are also necessary to implement an adequate migration (30,31).
Recent studies have revealed the importance of autophagy (an intracellular process that controls protein and organelle degradation and recycling) in the control of locomotion (32).The turnover of focal adhesions also regulates cell spreading and migration (33).Calcium ions (Ca 2+ ), which impact globally on almost every aspect of cellular life, play an important role in the control of directed movements (34).In this sense, the endoplasmic reticulum, a multifunctional signaling organelle which controls a wide range of cellular processes such as the entry and release of calcium ions, also participate in the regulation of cell locomotion (35,36).Cell polarity is required for an adequate directionality motion and there are a lot of molecular processes that have been implicated in the intrinsic polarity status of cells; in this regard, the centrosomes positioning serves as a steering device for the directional movement (37,38) and dynein together with other molecules regulates centrosomal orientation to establish and maintain cell polarity (39).
The Golgi apparatus (another important molecular processing center for modified proteins received from the endoplasmic reticulum) allows the remodeling of intracellular traffic processes toward the direction of movement; therefore, signals from the Golgi matrix play an important role in cell motility (40).The nucleus is very important for developing appropriate mechanical responses during cell migration, in fact this organelle behaves as a central mechanosensory structure, and its physical properties strongly connected to the cytoskeleton guarantee a proper cell migration (41,42).Recently, it has been described that structural chromatin organization also has a key role in the cellular migration process (43).In addition, many molecules and corresponding processes are involved in directional movement of cells as, for instance, focal adhesion proteins (talin, paxillin, vinculin, and others) (28), SCAR/WAVE proteins (44), actin-binding proteins (45), p21-activated kinases (a family of serine/threonine kinases) (46), TORC2/PKB pathway (47), mitogen-activated protein kinases (48), Arp2/3 complexes (49), WASP family proteins (50), Nck family of adaptor proteins (51), etc.In addition, numerous studies provide that integrative functional responses underlie in cell functionality which allow the emergence of systemic behaviors (52)(53)(54)(55)(56)(57).
All this evidence suggests that cell migration is not a mere metabolic-molecular process which can be regulated by any of its individually considered components.Most, if not all, fundamental cellular physiological processes appear to be involved in cellular locomotion, which is indicative of the emergence of a global functional phenomenon in the cell.However, confirming the systemic nature of cellular locomotion represents a scientific challenge of great difficulty.Such verification requires multidisciplinary approaches that combine complex experimental studies with advanced quantitative methods.
Here, we have addressed the key question: Is cell migration a highly coordinated and integrated emergent process at the global cellular level?To answer this question, we have designed a large quantitative study to analyze the systemic trajectories of 700 individual cells belonging to three different species: Amoeba proteus, Metamoeba leningradensis, and Amoeba borokensis.Such analysis has been performed under four different scenarios: in absence of stimuli, under chemotactic gradient (we have used an nFMLP peptide, which indicates to the amoebae the possible presence of food in their immediate environment), in an electric field (the electric membrane potential of cells enables predators like amoebas the detection of preys), and under complex external conditions such as simultaneous galvanotactic and chemotactic gradient stimuli.
To understand the forces driving the locomotion movement of the cell, all trajectories were analyzed using computational methods and advanced nonlinear physical-mathematical tools rooted in Statistical Physics (Statistical Mechanics).These quantitative studies focused on some essential characteristics of the systemic dynamics underlying locomotion movements.The results indicate that a very complex dynamic structure emerges in the migratory movements of all the cells analyzed.Such structure is mainly characterized by highly organized move-step sequences with very low entropy and high information, marked interdependence in the move steps with power-law autocorrelation decays, strong anomalous superdiffusion dynamics, persistence effects with trend-reinforcing behavior, and efficient movements to explore the extracellular medium.
This outstanding cellular dynamic structure is a consequence of the emergent systemic dynamics occurring in the cell.The locomotion movements seem to depend on a complex integrated selforganized system carefully regulated at global level, arising from the cooperative nonlinear interaction of most, if not all, cellular components.Such emergent systemic properties are not found specifically in any of the molecular parts, partial mechanisms, or individual processes of the cell.

Results
The migratory trajectories of 700 individual cells belonging to the three species, A. proteus, M. leningradensis, and A. borokensis, were recorded in four different scenarios: in absence of stimuli, under chemotactic gradient, in an electric field, and under simultaneous galvanotactic and chemotactic stimuli.Amoebae show robust movement in response to an electric field in a range between 300 and 600 mV/mm (galvanotaxis).Under such conditions, practically all amoebae migrate toward the cathode (58).Likewise, these cells also exhibit chemotactic movements.More specifically, the peptide nFMLP (N-formylmethionyl-leucyl-phenylalanine) secreted by bacteria may indicate that food might be in the near environment, provoking a strong chemotactic response (59).
All our experiments were performed in a specific setup consisting of two standard electrophoresis blocks (17.5 cm long), two agar bridges, a power supply, and in the middle of the experimental platform, a structure of standard glass slide and covers where the cells were located (see Fig. 1A-E and supplementary material).One electrophoresis block was directly plugged into a normal power supply and the other was connected to the first one through two agar bridges, thus preventing the direct contact of the anode and cathode with the medium (Chalkley's simplified medium (58)) where the cells were placed.Specifically, the amoebae were arranged in the center of the structure of standard glass slide and covers (experimental chamber) and their migratory displacements were monitored.The glass experimental structure enabled the generation of a laminar flux allowing the electric current to pass through, on one hand, and generating an nFMLP peptide gradient, on the other (supplementary material).
Prior to each experiment, all cells were starved for 24 h.The individual migratory movements of each cell were recorded over periods of 30 min using a digital camera attached to a stereo microscope.The experiments on flat 2D surfaces were always made with small groups of cells (no more than nine cells per replication).The following basic experimental information data (BEID) is provided for each scenario: "Nr" is the number of cells per replication, "Er" is the number of experimental replications, and "N" is the total number of cells.Finally, the recorded trajectories were analyzed in the form of time series using advanced nonlinear dynamic tools.

Cellular migratory movements without external stimulus
First, we recorded the locomotion trajectories of 153 individual cells belonging to the 3 species considered in a medium without any external influence (BEID: A. proteus: n = 50, Er = 7, Nr = 7-8; M. leningradensis: n = 51, Er = 7, Nr = 5-8; A. borokensis: n = 52, Er = 7, Nr = 6-8).In Fig. 2A, a representative example of these amoebae migratory movements in absence of stimuli is depicted (for clarity only 60 cells were randomly taken from the total).It can be observed that after 30 min, cells have explored practically all the directions of the experimentation chamber.To quantitatively analyze cell directionality, we calculated the displacement cosine for each trajectory, 153 cells in total (Fig. 2A′).Values close to −1 indicate a preference toward the left, while values close to 1 suggest a preference toward the right.Our analysis showed that values ranged between −1 and 1, with a median/IQR (interquartile range) of 0.05/1.40.Median/IQR values for each species were 0.42/1.24(A.proteus), −0.25/1.21(M.leningradensis), and 0.13/1.23 (A.borokensis).These results indicate that in absence of stimuli cells moved randomly without any defined guidance.

Cell migration under galvanotaxis conditions
The migratory trajectories of 147 cells belonging to the three species were recorded under an external controlled direct-current electric field of about 300-600 mV/mm.In Fig. 2B, a representative example of the migratory movements of 60 cells is depicted.They show an unequivocal systemic response consisting of the migration to the cathode which has been placed on the right side of the setup.The overall median/IQR value of the displacement cosines of all 147 cells (Fig. 2B′) was 0.99/0.07.This finding confirmed that a fundamental behavior characterized by an unequivocal directionality toward the cathode had emerged under these galvanotactic conditions.The median/IQR values for each species were 0.99/0.02(A.proteus), 0.98/0.10(M.leningradensis), and 0.99/0.08(A.borokensis).We compared the distributions of the values for the displacement cosines under galvanotaxis with the values obtained in the experiment without stimuli using the Wilcoxon rank-sum test.The results indicated that both behaviors were significantly different for the three species and that the galvanotactic cellular behavior is highly unlikely to be obtained by chance (P-values: 10 −9 , 10 −15 , and 10 −12 ; Z: −5.85, −7.79, and −6.85 for A. proteus, M. leningradensis, and A. borokensis, respectively).BEID: A. proteus: n = 49, Er = 7, Nr = 6-8; M. leningradensis: n = 48, Er = 7, Nr = 6-8; A. borokensis: n = 50, Er = 8, Nr = 3-9.

Cell locomotion under chemotaxis conditions
The migratory behavior of 166 cells belonging to the three species considered was recorded under conditions of chemotactic gradient.All the amoebae were exposed for 30 min to an nFMLP peptide gradient which had been placed on the left side of the setup.In Fig. 2C, a representative example of these migratory trajectories with 60 cells is depicted.Under these conditions 78.31% of all studied amoebae showed locomotion movements toward the attractant peptide.
The displacement angle cosines of the 166 individual trajectories ranged from −1 to 1, with a median/IQR value of −0.67/0.87.Median/IQR values for each species were −0.65/0.75(A.proteus), −0.77/0.68(M.leningradensis), and −0.51/1.17(A.borokensis), indicating that they exhibited a single fundamental behavior of movement toward the peptide (Fig. 2C′).The Wilcoxon rank-sum test showed significant differences between the cosine values obtained with and without chemotactic stimulus (P-values: 10 −7 , 0.00, and 0.02; Z: 5.08, 2.94, and 2.40 for A. proteus, M. leningradensis, and A. borokensis, respectively), and between the cosine values with chemotactic gradient and with the presence of an electric field (P-values: 10 −15 , 10 −17 , and 10 −14 ; Z: 8.00, 8.55, and 7.68 for A. proteus, M. leningradensis, and A. borokensis, respectively).This confirmed that the systemic locomotion behavior under the chemotactic gradient was completely different from both the absence of stimuli and the presence of an electric field.

Cellular displacement under simultaneous galvanotactic and chemotactic stimuli
Once the locomotion movements of the cells were recorded under the three previous independent experimental scenarios (without stimuli, under galvanotaxis, and under chemotaxis), we studied the trajectories of 234 cells under simultaneous galvanotactic and chemotactic stimuli.For such a purpose, the nFMLP peptide was arranged on the left of the setup (in the anode area) and the cathode was placed on the right.In Fig. 2D, a representative example of these locomotion movements (60 cells in total) is depicted.Under these complex external conditions, the results showed that 42% of the amoebae migrated toward the cathode while the remaining 58% moved toward the peptide (anode).
The displacement cosines of the 234 cells had an overall median/IQR value of −0.29/1.66.More specifically, the values for each species were (−0.32/1.59,median/IQR) for A. proteus, (−0.54/1.81,median/IQR) for M. leningradensis, and (−0.14/1.45,median/IQR) for A. borokensis.This analysis quantitatively verified that two main cellular migratory behaviors had emerged in the experiment, one toward the anode and another toward the cathode (Fig. 2D′).The statistical analysis (Wilcoxon rank-sum test) confirmed the presence of these two different behaviors for A. proteus (P-value = 10

Long-range interdependence in the move steps of cellular migratory displacements
An essential characteristic of systemic behavior in complex systems is the presence of dynamics with strong long-range correlations (60), "long-range" refers to decay in the autocorrelation function slower than exponential decay, where a single scale dominates the decay.One of the most recognized tools to analyze the presence of these correlations in time series (migratory trajectories here) is the "root mean square fluctuation" ("rmsf" analysis), a classical method in Statistical Mechanics based on the ideas raised by Gibbs (61) and Einstein (62).
Long-range interdependence can be detected by a power-law relation such that F(l) ∼ l α , where l is the number of steps.For uncorrelated data, the fluctuation exponent α is about 0.5, whereas α > 0.5 or α < 0.5 indicate, respectively, the presence of positive or negative long-range correlations (supplementary material).In Fig. 3A, an illustrative "rmsf" analysis for the locomotion movements of three representative cells belonging to each species considered under simultaneous chemotactic and galvanotactic stimuli (A.proteus and M. leningradensis) and under chemotactic conditions (A. borokensis) is depicted.
The results of "rmsf" analysis of the 700 experimental cell trajectories are shown in Fig. 3B (for more details, see Table S1).All the migratory trajectories exhibit long-range correlations in their cellular move-step migratory fluctuations.Specifically, we found that the scaling exponent α of the "rmsf" had a median/IQR value of 0.72/0.08 for A. proteus, 0.73/0.08 for M. leningradensis, and 0.71/ 0.08 for A. borokensis.The values of the "rmsf" analysis of the total experimental migratory trajectories analyzed ranged from 0.56 to 0.87, with a median/IQR of 0.72/0.08,whereas the values of the scaling exponent α of all shuffled trajectories ranged from 0.35 to 0.64, with a median/IQR value of 0.47/0.07(see Table S2 for more details).Moreover, a Wilcoxon test comparing measured exponents to those from shuffled trajectories revealed highly We also calculated the time duration of the correlations regime and found that all cells exhibited long-range correlations over periods ranging from 1.04 to 16.67 min with a median/IQR value of 9.38/7.29 min.These findings indicate strong dependences of past movements lasting approximately 1,125/875 (median/IQR) move steps (Fig. 3C and Table S3). A. proteus cells exhibited longrange correlations up to a median/IQR duration of 10.42/ 6.25 min, M. leningradensis cells showed 9.38/7.29 min, and A. borokensis cells 8.33/6.25 min, thus highlighting the influence of previous trajectory values on each cellular move step.These results show the presence of power-law autocorrelations decays in all migration trajectories.

Strong anomalous migratory dynamics in cellular locomotion
Another characteristic of the migratory movement of cells is their strong anomalous dynamics.This property is directly related to anomalous superdiffusion, a complex process with a high nonlinear relationship to time which also corresponds to efficient systemic directional trajectories (63,64).
One of the best methods to determine such dynamic property is the mean square displacement (MSD), a method proposed by Einstein (65) and later by von Smoluchowski (66).This Statistical Mechanics tool allows to quantify the amount of space explored by the amoebae during their locomotion.According to this procedure (see supplementary material), the anomalous diffusion exponent β is commonly used to refer to whether normal (Brownian, β = 1) or anomalous diffusion (β ≠ 1) is observed.The dynamics of subdiffusion and superdiffusion correspond to 0 < β < 1 and β > 1, respectively.
In Fig. 4A, we depicted an MSD analysis for the locomotion movements of three representative cells belonging to each species under absence of stimuli.The results of MSD analysis of the 700 experimental cells (see Fig. 4B and C and Table S4) show that practically all trajectories exhibit strong anomalous migratory dynamics.For experimental trajectories, the variable β, which characterizes the behavior of the diffusion process, had a median/IQR value of     and A. borokensis) in the four experimental scenarios ("Scenario 1" absence of stimuli, "Scenario 2" presence of an electric field, "Scenario 3" presence of a chemotactic peptide gradient, and "Scenario 4" simultaneous galvanotactic and chemotactic stimuli).A′) the percentage of cells moving in any specific sector of the experimental chamber is represented as a polar histogram divided into eight areas of π/4 angle amplitude each.B′-D′) histograms of the displacement cosines from panels B-D), respectively.E) Digitized cell trajectory with two inserts highlighting displacements regions of interest."N" is the total number of cells; "t" is the experiment duration; "p" is the chemotactic peptide (nFMLP); "+" represents the anode; "−" represents the cathode; "Sc1" Scenario One (absence of stimuli).Both the x-and y-axes show the distance in mm, and the initial location of each cell has been placed at the center of the diagram.

Complexity and information in cellular migration
To assess the information content within locomotion trajectories, we implemented the approximate entropy (ApEn), a robust approximation of the Kolmogorov-Sinai (K-S) entropy (67,68), providing insight into the complex migratory behavior that emerges from the cellular system.In Fig. 5A (Tables S6 and S7), the results of the ApEn estimation for the 700 cellular trajectories are shown.The heatmaps display the approximate K-S entropy for all experimental (upper row) and shuffled trajectories (bottom row) from each species, calculated for 72 different time windows (intervals) of increasing length  ) for all species and experimental conditions.C) Estimated distribution, median, and mean MSD β exponent values from experimental trajectories are illustrated using violin plots."Sc1" Scenario One, absence of stimuli; "Sc2" Scenario Two, presence of an electric field; "Sc3" Scenario Three, presence of a chemotactic peptide gradient; and "Sc4" Scenario Four, simultaneous galvanotactic and chemotactic stimuli.(interval duration was increased by 25 s at every iteration).
Intervals present ApEn values that vary from 10 −4 to 0.52 for experimental trajectories and from 0.14 to 2.13 for shuffled trajectories.These findings allow to observe that practically all the experimental series exhibit extremely low entropy.In Fig. 5B (Table S6), the decrease in entropy that occurs from SC1 to SC2-SC4 denotes how in absence of stimuli (SC1) cells maximize entropy, but when there is a more defined directionality cellular trajectories become more focused (less entropic displacements).Specifically, we found that the ApEn values of experimental trajectories exhibited a narrow range of low values displaying a median/IQR of 0.00/0.00for A. proteus, 0.00/0.00for M. leningradensis, and 0.00/0.00for A. borokensis.The ApEn analysis for all experimental migratory trajectories obtained under the four scenarios showed a median/IQR ApEn value of 0.00/0.00with values ranging from 10 −4 to 0.02.The ApEn values for shuffled trajectories displayed a range of very high values (from 1.25 to 2.13, median/ IQR equal to 1.97/0.16)relative to the values for experimental trajectories (Table S7).
The whole analysis confirms the presence of a complex structure characterized by high information in the move-step sequences in the migration trajectories of all cells.Furthermore, the statistical analysis revealed that this complex dynamic structure observed in the move-step trajectories was highly unlikely to occur by chance, as indicated by a P-value ≅ 0 and Z = −32.39,results of a Wilcoxon test comparing the respective ApEn value distributions of experimental and shuffled trajectories.

Persistence in cellular migratory movements
Persistence is another main characteristic of the systemic cellular migratory movements in unicellular organisms (69,70).The detrended fluctuation analysis (DFA) (see supplementary material) is a well-known technique for measuring persistent effects in physiological time series.
For a given observation scale ℓ, DFA calculates the function F(ℓ) to quantify the fluctuations of the time series around the local trend.If the time series displays scaling properties, then F(ℓ) ∼ ℓ γ asymptotically, where γ represents the scaling exponent.This exponent is commonly estimated as the slope of a linear fit in the log(F(n)) vs. log(ℓ) plot.Thus, γ serves as a measure persistence and helps to characterize the underlying dynamical system.Specifically, values close to 0.5 indicate the absence of long-range correlations, while when 1.5 < γ < 2, the process exhibits positive long-range persistence (71) (Fig. 6A).Through the application of this quantitative method, we identified the presence of long-range persistence in all experimental trajectories (Table S8), with a γ overall median/IQR value of 1.78/0.11.Specifically, the median/ IQR DFA scaling parameter γ was found to be 1.80/0.09for A. proteus, 1.79/0.14for M. leningradensis, and 1.78/0.10for A. borokensis (Fig. 6B and Table S8), thus indicating that all the move-step trajectories exhibit "trend-reinforcing behavior" (significant persistence).
In order to assess the reliability of the DFA analysis, we conducted a random shuffling procedure on 700 time series.The results demonstrated that the strong correlation values observed in the experimental migration series vanished after shuffling (refer to Fig. 6B and C and Table S9 for more information), with γ overall median/IQR of 0.48/0.14.This finding confirms that the complex locomotion structure, characterized by well-organized move-step sequences and persistent dynamics observed in the migration trajectories of the three cell groups, is not attributable to a random chance (P-value ≅ 0, Z = 32.39).

Kinematic properties in cellular locomotion trajectories
To quantify some kinematic properties of the cell migration trajectories, we studied the Intensity of the response (IR), the directionality ratio (DR), and the average speed (AS) of amoebae (Fig. 7A-C).
The IR is associated to the space explored by the cell, and in particular, we quantified the module of the trajectories to represent the strength of the response.In this case, the median/ IQR IR was 5.32/3.5 for A. proteus, 4.96/5.05for M. leningradensis, and 3.58/3.07for A. borokensis.Next, we studied the DR, which quantifies the trajectory straightness, ranging between 0 (for fully curved trajectories) and 1 (for fully straight trajectories), by considering the start and end point of the trajectory.The values ranged between 0.052 and 0.88 (median/IQR 0.53/0.31)for A. proteus, 0.04 and 0.87 (median/IQR 0.51/0.32)for M. leningradensis, and 0.01 and 0.95 for A. borokensis (median/IQR 0.47/0.30).
As it can be observed from the P-values derived from Kruskal-Wallis analyses there is a remarkable variability regarding kinetic properties, both between species (for example, the P-values comparing the IR, DR, and AS in Scenario 4 were 10 −13 , 10 −4 , and 10 −16 , respectively) and between scenarios (for example, the P-values of A. proteus for IR, DR, and AS compared among all four scenarios were 10 −6 , 10 −9 , and 10 −10 , respectively); for more information, see Fig. 7 and Table S10.
Figure 7D-G shows a clustering analysis performed on all kinematic properties considered.Each cluster was characterized by the proportion of cell types and experimental condition present in each group.The performance of the obtained clustering solution was assessed using the Silhouette coefficient, which estimates all the differences between intracluster points minus the distances between intercluster points.A higher Silhouette index indicates a model with better defined clusters.The implementation was achieved using the silhouette score implemented in Scikit-Learn.The three clusters identified in Fig. 7 yielded a Silhouette coefficient of 0.37.Similarly, the four clusters identified provided a Silhouette coefficient of 0.37.Another alternative clustering analysis by means of a hierarchical agglomerative method (Fig. S3) showed that the distinction between cell types or experimental conditions remains unchanged when varying the clustering strategy, which indicates the robustness of the findings.
The high variability of the statistics of the kinematic parameters, combined with the cluster analysis, indicates a high level of heterogeneity in any of the three metrics used.This heterogeneity exists both among cell types and among different experimental conditions and suggests that the behavior is individual in each one and hence they cannot be categorized into groups.

Dynamic structure in migratory movements
Finally, we have represented all the main metrics considered in our study, such as RMSF Alpha, RMSF correlation time (measured in move steps or in minutes), DFA Gamma, MSD Beta, and approximate entropy by comparing them with those data obtained in the corresponding shuffling procedures (Fig. 8A-C).As evident in each panel, the metrics' shuffled and nonshuffled values could be distinctly grouped and differentiated.This suggests that the inherent systemic information structure was disrupted during the shuffling process.The trajectories of the experimentally observed cells are completely differentiated from the cells whose trajectories lost systemic properties.
In Fig. 8G-J, a clustering analysis of the main metrics was performed.The three clusters related to cell type identified in Fig. 8 yielded a Silhouette coefficient of 0.35.Similarly, the four clusters identified related to experimental conditions provided a Silhouette coefficient of 0.31.This unsupervised clustering analysis combined to the small variability of the results, appears to be quite homogeneous, regardless the cell type or scenario considered, since its quantitative aspects show no dependency on either cell type or experimental condition.Moreover, we utilized an alternative approach, specifically hierarchical agglomerative clustering (Fig. S4), and the resulted profile remained consistent when the different clustering methods were applied, highlighting the robustness of the results.These findings suggest the emergence of a highly intricate dynamic structure within the migratory patterns of all examined cell trajectories.Furthermore, this structure appears to be an inherent aspect of cell locomotion, irrespective of species or environmental conditions.Indeed, cluster analysis of all quantitative parameters revealed no reliance on either cell type or experimental context, suggesting the potential universality of this behavior.

Discussion
Cellular migration is a cornerstone issue in many essential physiological and pathological processes.Here, we have addressed the integrative systemic dynamics involved in the regulation of directional motility.To this end, we have studied the

Amoeba proteus Metamoeba leningradensis Amoeba borokensis
Median Mean Fig. 5. Complexity and information in cellular migration.A) Heatmaps for the approximate entropy values of all 700 experimental (upper row panels) and shuffled (bottom row panels) cell trajectories from each species (A.proteus, M. leningradensis, and A. borokensis).Each row in every panel corresponds to a single cell, while in the 72 columns the endpoint of the approximate entropy calculation is represented, increased in 25 s at every iteration.B) Violin plots illustrate the estimated distribution, mean, and median approximate entropy values for all experimental cell trajectories."Sc1" Scenario One, absence of stimuli; "Sc2" Scenario Two, presence of an electric field; "Sc3" Scenario Three, presence of a chemotactic peptide gradient; and "Sc4" Scenario Four, simultaneous galvanotactic and chemotactic stimuli.migratory displacements of 700 single cells, belonging to three different species (A.proteus, M. leningradensis, and A. borokensis), in four different scenarios: in absence of stimuli, under chemotactic gradient, in an electric field and under complex external conditions such as simultaneous galvanotactic and chemotactic gradient stimuli.The experimental trajectories, obtained on flat 2D surfaces, have been quantitatively studied using a multidisciplinary approach to understand how integrative systemic forces drive the locomotion movement of cells.First, we have analyzed the interdependence in the move steps using the "rmsf" method, a classical Statistical Mechanics approach based on the concepts developed by Gibbs (61) and Einstein (62), later developed and used to quantify biological processes (72)(73)(74).Our results demonstrate that each move step ahead at a given point is strongly influenced by its preceding displacements, indicating that strong dependences of past movements lasting approximately 1,137 move steps over periods averaging 9.5 min do exist.Practically, all the 700 unicellular organisms analyzed in the four experimental conditions exhibited nontrivial correlations in their directional trajectories, which represents a key characteristic of the systemic dynamic movements emerging in the cell system (69).
Anomalous behavior is another characteristic of the migratory movements that we have identified using the MSD method also proposed by Einstein (65).Migratory dynamics that do not result in a linear MSD can be considered as nontrivial.Specifically, the anomalous nature of cell migration can be detected by superdiffusion, a physical phenomenon detected in the trajectories of all the 700 cells analyzed.Likewise, the MSD is a proxy for the surface area explored by the cell over time and is a measure related to the overall migration efficiency.The cellular displacements analyzed correspond to efficient movements during the exploration of the extracellular medium (63,64,75).The strong manifestation of the anomalous nature of cell migration can be caused by temporal memory effects as a consequence of the correlations in the cellular move steps (75)(76)(77)(78).
We have also quantified the regularity and unpredictability of the fluctuations over the migratory displacements (67,68).The obtained results show high levels of information in all the analyzed trajectories.This finding, together with the previous ones, confirms the presence of a very complex structure in the migratory move-step sequences.Entropy is directly related to the complexity of the system dynamics, and the very low level of entropy in the directional movements indicates that the migration patterns are organized on a level of complexity that is above the individual components of the cell system.Such complex dynamic structure observed in the trajectories was highly unlikely to occur by chance (practically, P-value ≅ 0).
In addition, we have verified the presence of persistent effects in the cellular migratory movements.The results of the DFA fluctuation analysis (79) show that the scaling exponent γ displays a total average (± SD) of 1.75 ± 0.12, indicating that the move-step trajectories exhibit a trend-reinforcing memory, that is, if the directional movements in the past show an increase in a set of their move-step values it is very likely to be followed by an increasing trend in the future; and vice versa, a decreasing trend in the past, is likely to be continued in the future.In other words, the evolution of the cell system trajectories is strongly influenced by previous system movements over long periods of time.Persistent effects are a key concept, closely related to temporal correlations, widely developed in Physics with a robust and formal Mathematical construction.Therefore, the results we have obtained with DFA also validate the presence of strong correlations in the locomotion movements.It is necessary to note that temporal correlations and regimes of anomalous diffusion have also been observed in the analysis of cellular trajectories previously (80)(81)(82)(83).
The quantitative studies carried out here unequivocally show that a very complex dynamic structure emerges in the migratory movements of all the analyzed cell systems.Such structure is characterized by highly organized move-step sequences with very low level of entropy and high information, nontrivial temporal interdependence in the move steps with power-law autocorrelation decays, strong anomalous superdiffusion dynamics, persistent effects with trend-reinforcing behavior, and efficient movements to explore the extracellular medium.The outstanding detected dynamic structure underlies all the migration trajectories of 700 cells of three different species analyzed under the four experimental scenarios.On the other hand, the results of the two types of clustering analysis performed suggested the potential universality of this complex systemic structure in the cellular locomotion movements.
These essential characteristics of the locomotion movements are a consequence of the self-organized dynamics intrinsic to all unicellular organisms.Cells are sophisticated systems conformed by the mutual interactions of millions of molecules and hundreds of thousands of macromolecular subcellular structures (57).They are open systems that operate far from the thermodynamic equilibrium and exchange energy-matter with the environment (84)(85)(86).Under these conditions, nonlinear enzymatic interactions G H J Fig. 8. Systemic metrics and clustering analysis of cellular migratory displacements.A-F) Analysis of the 700 experimental cell trajectories with the main metrics considered in our work such as RMSF alpha, DFA gamma, MSD beta, and approximate entropy is compared with corresponding shuffled data in which the systemic informational structure was completely lost.G-J) Clustering analysis of cellular motion experiments does not distinguish between cell types and experimental conditions when using nonlinear advanced movement metrics.Similar to Fig. 7D-G but defining each experiment using the 5D main vectors of metrics.k-Means clustering was performed using the first three principal components (PC1, PC2, and PC3), which accounted for 92.54% of the total variance.The interpretation of the different panels is the same as in Fig. 7, but the variables used for clustering are now different.
and irreversible metabolic processes allow the cell system to become spatially and temporally self-organized (56,57,(87)(88)(89)(90).If cells reach the equilibrium, their sophisticated dynamic functionality and molecular order disappears and they die.Briefly described, the essential energy-matter flow generates a negative entropy variation inside the cell which corresponds to an emergent positive increment in the information of the system (91).Such information increases the complexity, producing collective functional patterns, highly ordered macrostructures, and complex self-organized behaviors as for instance molecularmetabolic rhythms and spatial traveling waves (86,92).
These emergent nonequilibrium molecular dynamics supported by permanent energy dissipation (continuously exporting entropy to the external medium) are known as self-organized dissipative structures (93,94).The principles of self-organization through energy dissipation were conceived and developed by the Nobel Prize Laureate in Chemistry Ilya Prigogine (60).
Intensive studies over the last six decades have demonstrated that cells are very complex self-organized dissipative systems (57,85,90,91,(95)(96)(97) in which integrated processes and systemic properties at different levels of organization and complexity do appear.At a basic level, dissipative molecular behaviors emerge, for example, in shaping actin polymerization waves involved in the cytoskeleton activities during cell migration (98,99), in selforganized oscillations in actin networks (100), in myosin dynamics (101), in microtubular behavior (102), and in intracellular calcium rhythms (103).At the highest level, complex systemic properties such as directional mobility, integral growth, reproduction, sensitivity to the external medium, adaptive responses, and evolution do occur (57).To note, these strong emergent properties cannot be found in their individual molecular components or in their single molecular-metabolic processes (90).
Systemic dynamics are emergent integrative processes in all unicellular organisms, and such behaviors are a consequence of the self-organization of the biochemical system as a whole (104).From collective metabolic-molecular constituents, all of them interacting nonlinearly with each other, emerge basic coherent selforganized structures and functional ordered patterns which originate a cell system that increases at different levels its structural and functional complexity driven by energy dissipation and molecular information processing (57,88,89).A critical attribute of these dissipative self-organized systems is the interacting dynamics exhibiting long-range correlations (105)(106)(107).
In brief: I. We have presented a relevant number of first-rate biologic molecular experiments (Introduction section), carried out by independent research groups, which show that practically all fundamental main physiological processes of the cell are involved in cell migration.These studies provide clear experimental evidence that integrative functional responses underlie the cellular locomotion movements.Such conclusive set of current specific experiments, unequivocally prove that the main functional physiological structures of the cell are involved in the regulation of cell migration.II.Our quantitative studies show that a very complex dynamic structure emerges in the migratory movements of all the analyzed cells.Such structure is characterized by highly organized move-step sequences with very low level of entropy and high information, marked temporal interdependence in the move steps with power-law autocorrelation decays, strong anomalous superdiffusion dynamics, persistent effects with trend-reinforcing behavior, and efficient movements to explore the extracellular medium.Such characteristics correspond to critically self-organized systems.The locomotion trajectories change continuously, since they exhibit random magnitudes that vary over time, but these stochastic movements shape a dynamic structure whose defining characteristics are preserved in all the conditions analyzed.This movement structure corresponds to complex behavior belonging to a self-organized cell system, in which the emergent integrative dynamics drive the locomotion movement of cells.III.We have analyzed the migratory displacements of a large number of cells (700 single cells in total), belonging to three different species (one of them very evolutionarily separated from the other two), in four different scenarios: in absence of stimuli, under chemotactic gradient, in an electric field, and under complex external conditions such as simultaneous galvanotactic and chemotactic gradient stimuli (which act as opposite attractive poles in our experiments).The similar migratory characteristics observed may be extendable to other unicellular organisms and therefore these characteristics are possibly universal across cell types.The results of the two types of clustering analysis here performed also suggested the potential universality of this complex emergent structure in the cellular locomotion movements of cells.IV.The first-rate biologic molecular experiments here presented (Introduction section) and the dynamical studies performed suggest that cell migration is an emergent systemic property, governed by the functional integration of most if not all the cellular components.Cellular locomotion seems to be regulated by complex integrated physiological processes, carefully regulated at a systemic level, which depends on the cooperative nonlinear interaction of most, if not all, cellular components.Therefore, the properties responsible of this locomotor behavior are not found specifically in any of their singular molecular parts, partial mechanisms, or individual processes in the cell.This fact does not invalidate the importance of studying the influence of individual metabolicmolecular pathways on cell migration.
Cell migration is a central issue in many human physiological and pathological processes.We consider that new researches combining quantitative migratory systemic dynamics with molecular studies are crucial for the development of nextgeneration, efficient cellular therapies for migration disorders.

Materials and methods
Details of all methods are described in the supplementary material.

Fig. 1 .
Fig. 1.Experimental setup layout.A and D) Top and lateral views of the experimental setup.1: anode; 2: cathode; 3: agar + KCl bridges; 4: chemotactic peptide; 5: electrode used to probe the electric field; 6: stirrer used to properly mix the peptide; 7: experimental glass chamber, the arrow signals the trajectory and direction of the laminar flow.B) Top view of the glass pieces that compose the experimental glass chamber.8: 75 × 25 mm standard glass slide; 9: longitudinal trimmed glasses; 10: sliding lateral glasses; 11: central piece of glass underneath which the cells are placed.C) Top view of the experimental glass chamber.E) Axial section of the experimental chamber.12: flow sectional area.The experimental chamber can be opened and closed by longitudinally displacing #10, allowing to place or remove cells when open and establishing a laminar flow of medium through #12 when closed (see supplementary material for further details).F) Close-up of the experimental setup before experimentation, devoid of medium and cells.G) All elements sketched and described in A-E are shown in real conditions.

Fig. 2 .
Fig. 2.Representative migration trajectories of the three species under four experimental scenarios.A-D) Migration behavior of the three species (A.proteus, M. leningradensis, and A. borokensis) in the four experimental scenarios ("Scenario 1" absence of stimuli, "Scenario 2" presence of an electric field, "Scenario 3" presence of a chemotactic peptide gradient, and "Scenario 4" simultaneous galvanotactic and chemotactic stimuli).A′) the percentage of cells moving in any specific sector of the experimental chamber is represented as a polar histogram divided into eight areas of π/4 angle amplitude each.B′-D′) histograms of the displacement cosines from panels B-D), respectively.E) Digitized cell trajectory with two inserts highlighting displacements regions of interest."N" is the total number of cells; "t" is the experiment duration; "p" is the chemotactic peptide (nFMLP); "+" represents the anode; "−" represents the cathode; "Sc1" Scenario One (absence of stimuli).Both the x-and y-axes show the distance in mm, and the initial location of each cell has been placed at the center of the diagram.

Fig. 3 .
Fig. 3.Long-range interdependence in the move steps of cellular migratory displacements.A) Log-log plot of RMSF F vs. l step for a representative cell of each species.The slope was α = 0.84 for A. proteus, α = 0.85 for M. leningradensis, and α = 0.84 for A. borokensis, indicating the presence of strong long-range interdependence in the move steps of all of them.B) Diagram representing the values (and the overall average ± SD) of all the scaling exponents α from cells belonging to each species (A.proteus, M. leningradensis, and A. borokensis) under each experimental scenario (Sc1-Sc4).C) Violin plots showing the estimated distribution, median and average memory persistence values from cellular trajectories."Sc1" Scenario One, absence of stimuli; "Sc2" Scenario Two, presence of an electric field; "Sc3" Scenario Three, presence of a chemotactic peptide gradient; and "Sc4" Scenario Four, simultaneous galvanotactic and chemotactic stimuli.

Fig. 4 .
Fig. 4.Strong anomalous migratory dynamics in cellular locomotion.A) Graphics showing the value of the exponent β by fitting log-log plots of MSD as a function of the time interval τ, for eight prototypic cells of each species (A.proteus, M. leningradensis, and A. borokensis).β = 1 indicates normal diffusion, while β = 2 indicates ballistic diffusion.The gray region defines the area of superdiffusion, which is a complex process with a high nonlinear relationship to time, within which all the experimental values fall.B) Diagram representing all values of the β exponents (and the overall average ± SD) for all cells of the three species in each experimental scenario (Sc1-Sc4) experimental values in red and shuffled values in blue.The shuffling step extinguished the long-term correlation structure, causing the sharp division between the experimental and shuffled value distributions (P ≅ 0, Z = −32.39)for all species and experimental conditions.C) Estimated distribution, median, and mean MSD β exponent values from experimental trajectories are illustrated using violin plots."Sc1" Scenario One, absence of stimuli; "Sc2" Scenario Two, presence of an electric field; "Sc3" Scenario Three, presence of a chemotactic peptide gradient; and "Sc4" Scenario Four, simultaneous galvanotactic and chemotactic stimuli.
De la Fuente et al. | 7

Fig. 6 .
Fig. 6.Long-term memory effects in cellular migratory movements.A) Log-log plot of the detrended fluctuation parameter F(n) vs. window size n for a prototype cell from each species.The scaling exponent γ was γ = 1.90 for the A. proteus cell, γ = 1.90 for the M. leningradensis cell, and γ = 1.86 for the A. borokensis cell, indicating long-range memory effects in all the species.The shuffling procedure removed the memory information contained in the original trajectories, causing γ values to drop to γ = 0.51 for the representative A. proteus cell, γ = 0.56 for the M. leningradensis cell, and γ = 0.47 for the A. borokensis cell.B) Diagram displaying all the values of the scaling exponent γ (and the overall average ± SD) in all cells, separately for each species (A.proteus, M. leningradensis, and A. borokensis) and experimental scenario (Sc1-Sc4).Values of the scaling exponent γ belonging to shuffled time series are depicted in blue, while exponents corresponding to experimental trajectories are depicted in red.C) Violin plots illustrate the estimated distribution, mean, and median scaling exponent γ values for all experimental (upper row) and shuffled (bottom row) cell trajectories."Sc1" Scenario One, absence of stimuli; "Sc2" Scenario Two, presence of an electric field; "Sc3" Scenario Three, presence of a chemotactic peptide gradient; and "Sc4" Scenario Four, simultaneous galvanotactic and chemotactic stimuli.

Fig. 7 .
Fig. 7. Kinematic properties and clustering analysis of cellular locomotion trajectories.A-C) Group boxplots of the distributions of all experimental values for three main cytokinetic metrics, revealing basic data about the cell migration characteristics A) IR, B) DR, and C) AS of the three species used in this experiment (A.proteus, M. leningradensis, and A. borokensis) under all four experimental scenarios (Sc1, Sc2, Sc3, and Sc4).D-G) Unsupervised clustering was performed using k-means on experiments defined by 3D vectors of three metrics.D) Three clusters are depicted, with each dot representing a different cellular motion experiment (combining three cell types and four experimental scenarios).E) Similar to D but with four clusters.Note that the quantitative results of the cluster analysis in D are not the same as in E. F) Characterization of the three clusters (D) in terms of cell types and experimental conditions.G) Similar to F) but with four clusters (D).