Anomalous dynamics of cell migration
 *Institut für Physiologie, Medizinische Fakultät Carl Gustav Carus, Fetscherstrasse 74, D01307 Dresden, Germany;
 ^{‡}School of Mathematical Sciences, Queen Mary, University of London, Mile End Road, London E1 4NS, United Kingdom;
 ^{§}MaxPlanckInstitut für Plasmaphysik, EURATOM Association, D85748 Garching, Germany; and
 ^{¶}Institut für Physiologie II, RobertKochStrasse 27b, D48149 Münster, Germany
See allHide authors and affiliations

Edited by Charles F. Stevens, Salk Institute for Biological Studies, La Jolla, CA, and approved November 28, 2007 (received for review August 11, 2007)
Abstract
Cell movement—for example, during embryogenesis or tumor metastasis—is a complex dynamical process resulting from an intricate interplay of multiple components of the cellular migration machinery. At first sight, the paths of migrating cells resemble those of thermally driven Brownian particles. However, cell migration is an active biological process putting a characterization in terms of normal Brownian motion into question. By analyzing the trajectories of wildtype and mutated epithelial (transformed Madin–Darby canine kidney) cells, we show experimentally that anomalous dynamics characterizes cell migration. A superdiffusive increase of the mean squared displacement, nonGaussian spatial probability distributions, and powerlaw decays of the velocity autocorrelations is the basis for this interpretation. Almost all results can be explained with a fractional Klein–Kramers equation allowing the quantitative classification of cell migration by a few parameters. Thereby, it discloses the influence and relative importance of individual components of the cellular migration apparatus to the behavior of the cell as a whole.
Nearly all cells in the human body are mobile at a given time during their life cycle. Embryogenesis, woundhealing, immune defense, and the formation of tumor metastases are well known phenomena that rely on cell migration. Extensive experimental work revealed a precise spatial and temporal coordination of multiple components of the cellular migration machinery such as the actin cytoskeleton, cellsubstrate and cell–cell interactions, and the activity of ion channels and transporters (1–4). These findings are the basis for detailed molecular models representing different microscopic aspects of the process of cell migration such as the protrusion of the leading edge of the lamellipodium, or actin dynamics (5). Mathematical continuum models, in contrast, focus on collective properties of the entire cell to explain requirements for the onset of motion and some typical features of cell motility (6). These models are usually limited to small spatiotemporal scales. Therefore, they provide little information about how the integration of protrusion of the lamellipodium, retraction of the rear part, and force transduction onto the extracellular matrix lead to the sustained longterm movement of the entire cell. This process is characterized by alternating phases of directed migration, changes of direction, and polarization. The coordinated interaction of these phases suggests the existence of intermittency (7) and of strong spatiotemporal correlations. It is therefore an important question whether the longterm movement of the entire cell can still be understood as a simple diffusive behavior like usual Brownian motion (8, 9) or whether more advanced concepts of dynamic modeling have to be applied (10, 11).
Results and Discussion
We performed migration experiments and analyzed the trajectories of two migrating transformed renal epithelial Madin–Darby canine kidney (MDCKF) cell strains: wildtype (NHE^{+}) and NHEdeficient (NHE^{−}) cells (12). The cells were observed for up to 1,000 min. Fig. 1 a depicts the contours and the path of a migrating MDCKF NHE^{+} cell monitored for 480 min. At first sight, the cell's trajectory resembles those of normal Brownian particles. Brownian motion in terms of the Ornstein–Uhlenbeck process (13, 14) is characterized by a mean squared displacement, msd (for definition, see Eq. 1 ), proportional to ∼t ^{2} at short times corresponding to ballistic motion and ∼t for long time intervals designating normal diffusion. Our experimental data show that both types of MDCKF cells behave differently. Consistent with earlier observations, MDCKF NHE^{−} cells move less efficiently than NHE^{+} cells (12, 15, 16), resulting in a reduced msd for all times. As displayed in Fig. 1 b, their msd exhibits a crossover between three different dynamical regimes. For short times (<4 min; phase I), the increase of the msd differs from a ballistic t ^{2} scaling. The logarithmic derivative of the msd (Eq. 2 ) shown in Fig. 1 c characterizes this first region with an exponent β(t) below ≈1.8. In the subsequent intermediate phase II (up to ≈20 min), the msd reaches its strongest increase with a maximum exponent of β(t) ≈ 1.8. When the cell has approximately moved beyond a squared distance larger than its mean squared radius, r (indicated by arrows in Fig. 1 b), the exponent β(t) of the msd gradually decreases below 1.5.
We next extracted the probability that the cells reach a given position x at time t from the experimental data. This corresponds to the temporal development of the spatial probability distribution function p(x, t) delivering information beyond the msd. Fig. 2 a and b reveals the existence of nonGaussian p(x, t) distributions for different points in time. The transition from a peaked distribution at short times t = 1 min to rather broad distributions at long times t = 120 min and t = 480 min in Fig. 2 a and b suggests again the existence of distinct dynamical processes acting on different time scales. The shape of the distributions can be quantified by calculating the kurtosis according to Eq. 5 , which is displayed as a function of time in Fig. 2 c. The kurtosis rapidly decays from values of ≈7–9 to reach a constant value of ≈2.3 for both cell types in the long time limit. Such a behavior implies a transition of p(x, t) from a peaked to a flat form also visible from Fig. 2 a and b. The timedependent deviation of the kurtosis from the value of 3, which would correspond to a Gaussian distribution (as, e.g., for the Ornstein–Uhlenbeck process), is another strong manifestation of the anomalous nature of cell migration.
To gain further insight into the origin of the anomalous dynamics, we calculated the velocity autocorrelation function v _{ac}(t) (defined in Eq. 4 ) that characterizes the correlation of the velocity (Eq. 3 ) at time t with its value at time t = 0 (Fig. 3). As for the msd in Fig. 1 b, the doublelogarithmic plot displays transitions between three regimes characterized by different time scales for both cell types. After a pronounced dip at short times, the autocorrelation function shows a gradual transition during intermediate times to a powerlawlike decay at long times. In contrast, the Ornstein–Uhlenbeck realization of Brownian motion would imply a purely exponential decay of the velocity autocorrelation function for all times. The nontrivial behavior of the velocity autocorrelation function stresses the existence of longrange correlations in time of the underlying dynamical process of cell migration. We would like to emphasize at this point that no chemotactic or other physical gradients were imposed on the cells examined in our study. Thus, our analysis shows that the socalled “random migration” does actually not proceed as randomly as one might expect.
To generalize the interpretation of our data, we sought an integrative mathematical model. The experimental results posed extensive constraints on the choice of such a model. A full theory should generate a powerlaw behavior of the msd including transitions between different time scales, nonGaussian probability distributions, and a velocity autocorrelation function with a powerlaw decay for long times. Conventional Langevin, Fokker–Planck, and Klein–Kramers equations (17) provide descriptions of ordinary Brownian motion [e.g., the Ornstein–Uhlenbeck process (13, 14)] that do not match these constraints. The anomalous nature of cell migration demands for the inclusion of temporal memory in the above equations that can be achieved by introducing fractional derivatives (18–21). We therefore modeled our msd data with the fractional Klein–Kramers equation in Eq. 6 (18), which includes the Ornstein–Uhlenbeck process as special case (α = 1). The inclusion of an uncorrelated noise term delivers the formula for the msd as given in Eq. 10 . The continuous blue and orange lines in Fig. 1 b represent the resulting fits with parameter values as given in Table 1 for MDCKF NHE^{+} and NHE^{−} cells. There is excellent agreement of data and model for all times. The shallow initial slope of the msd curve is due to the estimated noise level η ≈ 0.78 μm and 0.47 μm for NHE^{+} and NHE^{−} cells, respectively. The values of η are much larger than the measurement uncertainty of σ_{pos} ≈ 0.1 μm showing the influence of “biological noise” generated by lamellipodial activity. The time scale τ_{α} = (1/γ_{α})^{1/α} (Table 1 and Eqs. 6 – 8 ) characterizes the transition of the MittagLeffler function from stretched exponential to powerlaw behavior for t ≫ τ_{α}. The resulting time scales τ_{α} = 18.7 min for NHE^{+} and τ_{α} = 15.2 min for NHE^{−} are comparable with the times at which the cells cross their mean squared radii (458 μm^{2} and 211 μm^{2} for NHE^{+} and NHE^{−} cell, respectively). For larger times, the msd shows a transition to a power law of ∼t ^{2−} ^{α} (see Eq. 8 ) thus describing superdiffusion of ∼t ^{1.25} for NHE^{+} and ∼t ^{1.28} for NHE^{−} cells. The enhanced thermal velocity, v _{th} ^{2}, and the slightly reduced value of γ_{α} generate a diffusion coefficient for NHE^{+} cells that is twice as large as that for NHE^{−} cells. Thus, our model quantitatively confirms the importance of the Na^{+}/H^{+} exchanger for directed migration (15, 16). The logarithmic derivative of the model Eq. 10 in Fig. 1 c emphasizes the influence of the noise term for short times generating the deviation from a ballistic initial increase of the msd in agreement with the experimental data. Without the noise term (η = 0), β(t) would be 2 for short times. In addition, Fig. 1 c shows a continuous transition of β(t) to the estimated exponent 2 − α (Table 1) for long times.
To show that the predictions of the fractional Klein–Kramers equation correspond more closely with the experimental data than those from simpler dynamical models, we included a quantitative analysis of the Ornstein–Uhlenbeck process (13, 14, 17). The application of Eq. 11 to the experimental msd data delivers the Ornstein–Uhlenbeck parameters in Table 2. At first sight, the doublelogarithmic plot in Fig. 1 b shows only small differences between the fractional model and the conventional Ornstein–Uhlenbeck process. However, the logarithmic derivative shown in Fig. 1 c is clearly in favor of the fractional Klein–Kramers equation, especially in phases II and III. For longer times, the Ornstein–Uhlenbeck process shows a too fast decay of the exponent β(t) toward values near 1 corresponding to normal diffusion.
The probability distribution p(x, t) of the fractional Klein–Kramers equation is only known in the limit of long times (t ≫ τ_{α}) given by the solution of the corresponding fractional diffusion equation (19) in terms of Fox functions (Eq. 9 ). Fig. 2 a and b compares these model solutions (using the parameters of Table 1) with data for times t = 120 and 480 min for MDCKF NHE^{+} and NHE^{−} cells, respectively. There is a good agreement between data and model for longer times. The rather peaked solutions in Fig. 2 a and b for the shortest time t = 1 min cannot be explained with these functions. For comparison, we have also added the Gaussian probability distributions of the Ornstein–Uhlenbeck process (green lines in Fig. 2 a and b), which can explain neither the peaked short time nor the long time behavior. Fig. 2 c illustrates that the kurtosis of the fractional Klein–Kramers equation deviates for short times, too. However, for longer times, data and theoretical kurtosis converge in agreement with the observation of p(x, t) and the Fox functions in Fig. 2 a and b toward a value of ≈2.3. In contrast, the Ornstein–Uhlenbeck process has a constant value of 3 over the entire time range.
Finally, in Fig. 3 we compare the velocity autocorrelation function of the fractional Klein–Kramers equation with the cell data. On the doublelogarithmic plot, the MittagLeffler function given by Eq. 7 (using the same parameters of Table 1 from the msd fit) nicely interpolates between the experimental data values, again showing the crossover between stretched exponential and powerlaw behavior within the same time scales as the msd in Fig. 1 b. In contrast, the velocity autocorrelation function of the Ornstein–Uhlenbeck process (green lines in Fig. 3) decays too fast, ∼exp(−γ_{1} t) (see Eq. 7 ), thereby missing the longrange correlations of the fractional model. The obvious deviations of the first data point at t = 1 min for both cell types from the fractional Klein–Kramers model can be explained with the contribution of the uncorrelated “biological noise” to the velocity autocorrelation function as given by Eq. 12 . Using the values of η from Table 1, the differences between the model solutions and the data points at t = 1 min are estimated as 1.23 μm^{2}/min^{2} and 0.45 μm^{2}/min^{2} for MDCKF NHE^{+} and NHE^{−} cells, respectively. This agrees quite well with the observed differences in Fig. 3 and is also the case for the corrections at t = 0 min, which are not visible in the doublelogarithmic plots.
In summary, we have shown that a variety of anomalous dynamical properties characterizes the migration process of MDCKF cells. In all of these quantities, we observe a crossover between anomalous dynamics on different time scales that reminds one of intermittent behavior as claimed to be important for optimal search strategies of foraging animals (7). The fractional Klein–Kramers equation amended by an uncorrelated noise term models the msd and velocity autocorrelation function for all times as well as the long time dynamics of the probability distribution p(x, t). Thus, our approach offers a theoretical framework that allows the classification of the dynamics of cell migration with a few physical parameters that can be calculated from the cells' trajectories. These can be compared for different cell types and under different experimental conditions. We probed our model by comparing migration of wildtype and NHEdeficient MDCKF cells. The defect in directional migration induced by NHE deficiency (12, 15) can clearly be detected by analyzing the dynamics of MDCKF cell migration. However, it is remarkable that the general pattern of anomalous dynamics is not changed by NHE deficiency. This observation indicates that the dynamics of cell migration is organized on a level of complexity that is above that of individual components of the cellular migration machinery. However, our model may allow the identification of those components of the cellular migration apparatus that govern the longterm behavior of migrating cells.
Materials and Methods
Cell Culture.
Experiments were carried out on transformed Madin–Darby canine kidney (MDCKF) cells and on NHEdeficient MDCKF cells (12) referred to as NHE^{+} and NHE^{−} MDCKF cells. Cells were kept at 37°C in humidified air containing 5% CO_{2} and grown in bicarbonatebuffered minimal essential medium (MEM) (pH 7.4) with Earle's salts (Biochrom) supplemented with 10% FCS (Biochrom).
Migration Experiments.
Cells were seeded at low density (to avoid collisions between cells during the experiment) in tissue culture flasks (Falcon) 1–2 days before the experiments. They were incubated in Hepesbuffered (20 mmol/liter, pH 7.4) MEM supplemented with 10% FCS during the course of the experiments. The flasks were placed in a heating chamber (37°C) on the stage of an inverted microscope equipped with phasecontrast optics (Axiovert 40; Zeiss). Migration was monitored with a video camera (Hamamatsu) controlled by HiPic software (Hamamatsu). Images were taken in 1min time intervals during a time range of up to 1,000 min. N = 13 trajectories of each cell type were used for the analysis consisting of >10,000 data points for each group.
Data Analysis.
Image segmentation was performed with Amira software (Visage Imaging/Mercury Computer Systems). The outlines of individual cells at each time step were marked throughout the entire image stack and taken for all further processing. In addition, we assessed the accuracy of the segmentation by repeated segmentation. The normally distributed experimental uncertainty amounts to σ_{pos} ≈ 0.1 μm (data not shown). Quantitative data analysis and calculation of parameters were performed with programs developed by us. The x and y coordinates of the cell center (in micrometers) were determined as geometric means of equally weighted pixel positions within the cell outlines as function of time. The combined trajectories of a cell population allow the calculation of the mean squared displacement, msd, that describes the mean of the squared distances between a common starting point at time t _{0} and the actual positions of a cell population at time t, where 〈…〉 denotes a combined average over all starting times t _{0} and cell paths. The increase of the msd can be quantified by the logarithmic derivative leading to a timedependent increase msd(t) ∼ t ^{β} ^{(t)}.
The velocity of migrating cells in x and y direction [v _{x/y}(t) (μm/min)] was calculated from trajectories as the difference quotient of two cell positions at times t + Δ and t with the observation time interval of Δ = 1 min. These velocities were used to calculate the velocity autocorrelation function where 〈…〉 denotes an average as in Eq. 1 .
The position of the cells can also be used to calculate the probability p(x, y, t) of finding a cell at position (x, y) for time t. Because we did not find any correlations between x and y direction in the velocity autocorrelation function (data not shown), we reduced the discussion to p(x, t) given as average of x and y positions.
The kurtosis is defined as the ratio of moments by It can be interpreted as shape index of the probability distribution function p(x, t) and takes the value of 3 for Gaussian functions.
Fractional Klein–Kramers Equation.
The anomalous properties of cell dynamics were assessed with the fractional Klein–Kramers equation (FKK) for the probability distribution P(x, v, t) in position x, velocity v, and time t as proposed by Barkai and Silbey (18) but without external forces: γ_{α} denotes the damping term, k _{B} is the Boltzmann constant, T is the temperature, M is the mass of the particle, and α defines the order of the fractional time derivative of Riemann–Liouville type (see, for example, ref. 20). For α = 1, the above equation reduces to the ordinary Klein–Kramers equation (17). The fractional Klein–Kramers equation of Eq. 6 implies that the velocity correlation function is given by the socalled MittagLeffler function (22) E _{α}: In the case of the Ornstein–Uhlenbeck limit (α = 1), E _{α} reduces to an exponential decay ∼exp(−γ_{1} t). The mean squared displacement of the fractional Klein–Kramers equation is represented by the generalized MittagLeffler function (22) E _{α,β}: The mean thermal velocity v _{th} ^{2} = k _{B} T/M is related to the generalized diffusion coefficient D _{α} by the relation D _{α} = v _{th} ^{2}/γ_{α}. In the limit of large damping (γ_{α} → ∞), the fractional Klein–Kramers equation reduces to a fractional diffusion equation (FDE) (19, 20) where the spatial probability distribution functions p(x, t) are given by a Fox function H (19, 20): Eq. 8 can be extended to include uncorrelated noise of variance η^{2} generated by measurement errors (23) or by biological activity; e.g., by the fluctuating lamellipodium (t > 0) Eq. 10 reduces to the Ornstein–Uhlenbeck result for α = 1 with the noise term In a similar way, uncorrelated noise also influences the velocity autocorrelation function. The fluctuations η_{i} affect the velocity calculation in Eq. 3 at both positions x(t _{i}) and x(t _{i} + Δ) separated by the measurement interval Δ. Performing the average for the calculation of the velocity autocorrelation function eliminates all linear terms η_{i} while quadratic noise terms deliver a modification of the velocity autocorrelation function at times t = 0 and t = Δ: with the Kronecker delta δ(t, t′) = 1 for t = t′ and 0 elsewhere. All other times are unaffected, if the noise source is uncorrelated.
Bayesian Data Analysis.
The parameters of the FKK model in Eq. 6 and their uncertainties were estimated with Bayesian data analysis (24) (for a recent review, see ref. 25) applied to the corresponding twodimensional mean squared displacement of Eq. 10 (or Eq. 11 for the Ornstein–Uhlenbeck process). Bayesian data analysis offers a logically consistent link between data and models. It allows a reliable estimation of the model parameters taking the uncertainties of the experimental data into account. The expectation value of the nth moment of the parameters θ̂ = {α, γ_{α}, v _{th} ^{2}, η} in Eq. 10 is given by integrating over the posterior probability function for these parameters p(θ̂  data). The later is proportional to the product of likelihood p(data  θ̂) and prior function p(θ̂): Assuming normally distributed errors, the likelihood function p(dataθ̂) is given by where msd(θ̂, t _{i}) corresponds to Eq. 10 for the FKK model. A constant prior p(θ̂) was applied because of the lack of advance information about the values of the parameters. The fourdimensional integral in Eq. 13 was evaluated numerically by Markov chain Monte Carlo sampling. Thus, application of Eq. 13 with n = 1 delivers the expectation values of the four parameters α (1), γ_{α} (1/min^{α})_{,} v _{th} ^{2} (μm^{2}/min^{2}), and η (μm) based on the available experimental msd data at time points t _{i}. The error σ(t _{i}) was estimated as where N is the number of cell paths. The quotient of path length T _{path} (= 500 min) and actual time t _{i} estimates the number of more or less independent measurement intervals during the calculation of the combined expectation value in Eq. 1 [also see Qian et al. (26) for the discussion of statistical errors of the msd]. In addition, Eq. 13 allows for the estimation of the uncertainties of the parameter with n = 2 via In addition, we applied this formalism to simulated data of the Ornstein–Uhlenbeck process with a number of data that are comparable with the experiments. These simulations showed that the presented Bayesian data analysis delivers an agreement of estimated and actually used simulation parameters within the uncertainties (data not shown).
Acknowledgments
We thank Eric Lutz for helpful discussions about anomalous dynamics, Andreas Deussen and Hans Oberleithner for fruitful comments, and Sabine Mally for excellent technical assistance. R.K. thanks the MPIPKS Dresden for frequent hospitality. This work was supported by Deutsche Forschungsgemeinschaft Grants Schw 407/93 and Schw 407/101 (to A.S.) and by British EPSRC Grant EP/E00492X/1 (to R.K.).
Footnotes
 ^{†}To whom correspondence should be addressed. Email: peter.dieterich{at}tudresden.de

Author contributions: P.D., R.K., R.P., and A.S. designed research, performed research, analyzed data, and wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.
 © 2008 by The National Academy of Sciences of the USA
References

↵
 Jalali S ,
 del Pozo MA ,
 Chen K ,
 Miao H ,
 Li Y ,
 Schwartz MA ,
 Shyy JY ,
 Chien S
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵

↵
 Stokes CL ,
 Lauffenburger DA ,
 Williams SK
 ↵
 ↵

↵
 Schwab A ,
 Rossmann H ,
 Klein M ,
 Dieterich P ,
 Gassner B ,
 Neff C ,
 Stock C ,
 Seidler U
 ↵
 ↵

↵
 Denker SP ,
 Barber DL

↵
 Stock C ,
 Gassner B ,
 Hauck CR ,
 Arnold H ,
 Mally S ,
 Eble JA ,
 Dieterich P ,
 Schwab A

↵
 Risken H

↵
 Barkai E ,
 Silbey R
 ↵
 ↵
 ↵

↵
 Erdelyi A
 ↵

↵
 Jaynes ET ,
 Bretthorst GL
 ↵
 ↵