Motile parameters of cell migration in anisotropic environment derived by speed power spectrum fitting with double exponential decay

i

This article has been updated. A Correction to this article was published on .

i

The Original Article was published on

i

The Original Article was published on

i

This article has been updated

i

This article was retracted on

i

This article has been updated

Abstract

Cell migration through anisotropic microenvironment is critical to a wide variety of physiological and pathological processes. However, adequate analytical tools to derive motile parameters to characterize the anisotropic migration are lacking. Here, we proposed a method to obtain the four motile parameters of migration cells based on the anisotropic persistent random walk model which is described by two persistence times and two migration speeds at perpendicular directions. The key process is to calculate the velocity power spectra of cell migration along intrinsically perpendicular directions respectively, then to apply maximum likelihood estimation to derive the motile parameters from the power spectra fitting with double exponential decay. The simulation results show that the averaged persistence times and the corrected migration speeds can be good estimations for motile parameters of cell migration.

Graphical abstract

Keywords

cell migration / heterogeneity / power spectrum / random walk / Langevin equation

1 Introduction

Cell migration is a basic function that is essential for many physiological and pathological processes [1], including embryogenesis [2, 3], nervous development [4], wounding healing [5, 6], inflammation [7] and immunological responses [8]. Controlled by complex cellular signaling pathways [9-11], the onset of ill-regulated cell motility is involved in many human diseases, among which cancer is the most typical.

As an essential biological phenomenon, cell motility has also gained extensive attention of physicists and mathematicians for a long time and has been studied as a random walk process [12-14]. In the absence of symmetry-breaking gradients [15], the motility of bacteria and eukaryotic cells has long been described by isotropic Brownian movement [16, 17]. As an inspired model, the persistent random walk (PRW) has been widely used to describe the isotropic motility of cells in 2D microenvironment [18-20], which is governed by the following Langevin equation [21, 22]:

\begin{array}{}\text{(1)}& \frac{\text{d}v}{\text{d}t}=-\frac{1}{P}v+\frac{S}{\sqrt{P}}\cdot \tilde{W}\mathrm{,}\end{array}

where \tilde{W}\sim N(\mathrm{0,1}) is Gaussian white noise, v the velocity vector, P the persistence time of cell migration and S the migration speed of the Ornstein-Uhlenbeck process [23].

With Eq. (1), the 2D cell migration is fully characterized by the motile parameters of migration speed and persistence time [24]. At the single-cell level, the persistence time quantifies the memory intensity of cell to past velocities. The PRW model can produce several featured properties of cell migration in the isotropic environment, including a Gaussian distribution of velocities, an exponential decay of the velocity auto-correlation function (VACF), an isotropic velocity field, and a linear mean squared displacement on a long time scale.

However, different from the normal diffusion, the anomalous diffusion migration has been found in many different kinds of cells in experiment. For example, non-Gaussian diffusion is more prevalent for particles undergoing the thermal jiggling in soft materials [25]. Regulated by chemokine CXCL10, a generalized Lévy walk is a better choice for CD8+ T cells to find rare targets in a smaller time [26]. Predators have higher encounter rates when adopting Lévy-type foraging in natural-like prey fields compared with purely random landscapes [27]. The distribution of cell velocities following a simple exponential function is found in the long-term cell migration in low-density monolayer cultures [28]. Tsallis’ distribution of velocity has also been proposed for endodermal Hydra cells in cellular aggregates [29]. The double-exponentially decaying VACF is discovered in the analysis of motility of human keratinocytes and fibroblasts [30]. Recent experiments showed that the metastatic breast cancer cells invade collectively the 3D collagen matrix by exchanging leaders in the invading front [31, 32], and the motility of fibrosarcoma cells through 3D extracellular matrices does not follow a random walk [33, 34].

As a consequence, different kinetic models have been proposed to explain different kinds of migration behaviors of cells [35-38]. Some of them are simply modified from the PRW model [30, 33, 34, 39], while others are essentially different from the PRW model [28, 32, 35, 36, 40-42]. Among those modified PRW models, the anisotropic PRW (APRW) is representative to describe the cell migration in anisotropic collagen matrix [33]. The APRW model can describe some novel properties, including a better fit of distribution of displacements, a double-exponentially decaying VACF, an anisotropic angular velocity magnitude and an uneven distribution of angular displacements on long time scales.

To explain cell migration behaviors by modeling, another research interest is to apply different methods to characterize quantitatively the migration trajectories of cells. For the widely-used PRW process, it is necessary to derive from the trajectories the two basic motile parameters, i.e., persistence time P and migration speed S. These two motile parameters can be fitted out with the discussion of mean squared displacement (MSD) [39, 43-47] and VACF [48] of migration trajectories. However, it has been pointed out that the fitting of P and S with MSD or VACF may return unreliable errors [49-51]. Considering that the power spectrum of velocity is the Fourier Transform [52, 53] (FT) of VACF [54, 55], the motile parameters can be obtained with least squared fitting or maximum likelihood estimation with FT of migration trajectories [56].

Different from the isotropic migration characterized by only two motile parameters in PRW model, the APRW model for 2D anisotropic migration has four motile parameters, i.e., two persistence times P and two migration speeds S at two perpendicular axes. In Ref. [57], a method has been proposed to obtain the two persistence times by fitting VACF from APRW trajectories under the condition that the two migration speeds are always equal for cell migration. However, it remains unknown how to derive reliably all the four motile parameters with APRW trajectories in general conditions. In this paper, we propose a method to solve this problem. We first distinguish different types of APRW processes based on different speed-acceleration profiles and the degree of heterogeneity of motile parameters. Then, we develop a method to obtain the four motile parameters with the application of singular vector decomposition [58, 59], power spectrum analysis [56, 60] and maximum likelihood estimation [49, 61]. The simulation results indicate that our method is effective and general to obtain the motile parameters of cell migration which can be described by APRW process.

2 Model and method

2.1 Anisotropic persistent random walk model

In APRW model [33, 34], cell motility is assumed to display different persistence times and migration speeds along two perpendicular axes. We regard the perpendicular directions as x and y axes in computer simulations, respectively. Then the velocities of cell migration along x and y axes are governed by the following Langevin equations:

\begin{array}{}\text{(2)}& \frac{\text{d}{v}_x}{\text{d}t}=-\frac{1}{P_x}{v}_x+\frac{S_x}{\sqrt{P_x}}\cdot \tilde{W}\mathrm{,}\end{array}

\begin{array}{}\text{(3)}& \frac{\text{d}{v}_y}{\text{d}t}=-\frac{1}{P_y}{v}_y+\frac{S_y}{\sqrt{P_y}}\cdot \tilde{W}\mathrm{.}\end{array}

In which P_x and P_y are persistence times, S_x and S_y are migration speeds. For a given set of parameters of P_x, S_x, P_y and S_y, the displacements of cell location in each time step of \text{d}t are governed by the following equations:

\begin{array}{}\text{(4)}& \text{d}x(t\mathrm{,}\text{d}t)={\alpha }_x\cdot \text{d}x(t-\text{d}t\mathrm{,}\text{d}t)+F_x\cdot \tilde{W}\mathrm{,}\end{array}

\begin{array}{}\text{(5)}& \text{d}y(t\mathrm{,}\text{d}t)={\alpha }_y\cdot \text{d}y(t-\text{d}t\mathrm{,}\text{d}t)+F_y\cdot \tilde{W}\mathrm{.}\end{array}

Here, \text{d}x and \text{d}y are the displacements of cell location in x and y axes in time step of \text{d}t. The terms {\alpha }_x, F_x and {\alpha }_y, F_y describe the migration behaviors of cell in x and y axes, respectively, defined as follows:

\begin{array}{}\text{(6)}& {\alpha }_x=1-\frac{\text{d}t}{P_x}\mathrm{,}{\alpha }_y=1-\frac{\text{d}t}{P_y}\mathrm{,}\end{array}

\begin{array}{}\text{(7)}& F_x=\sqrt{\frac{S_x^2\text{d}{t}^3}{P_x}}\mathrm{,}{F}_y=\sqrt{\frac{S_y^2\text{d}{t}^3}{P_y}}\mathrm{.}\end{array}

When taking the limit for persistence time, cell migration ranges from straight motion (P\to \infty ) to random walk (P\to \text{d}t), so the cell movement is the superposition of straight motion and random walk. When the motile parameters in x and y axes are the same, the APRW model becomes the isotropic PRW.

In this work, we define the heterogeneity index (HI) on the basis of the heterogeneity of motile parameters, as follows:

\begin{array}{}\text{(8)}& {\text{Φ}}_P=\max \left\{\frac{P_x}{P_y}\mathrm{,}\frac{P_y}{P_x}\right\}\mathrm{,}{\text{Φ}}_S=\max \left\{\frac{S_x^2}{S_y^2}\mathrm{,}\frac{S_y^2}{S_x^2}\right\}\mathrm{,}\end{array}

\begin{array}{}\text{(9)}& {\text{Φ}}_D=\max \left\{\frac{P_x{S}_x^2}{P_y{S}_y^2}\mathrm{,}\frac{P_y{S}_y^2}{P_x{S}_x^2}\right\}\mathrm{.}\end{array}

Here, {\text{Φ}}_P and {\text{Φ}}_S indicate the heterogeneity degree of persistence time and migration speed, respectively, and {\text{Φ}}_D represents the anisotropy degree of migratory ability. Note that all the values of {\text{Φ}}_P, {\text{Φ}}_S and {\text{Φ}}_D are greater than or equal to 1. Furthermore, we define it as positive anisotropy if both the persistence time and migration speed in one axis are larger than those in another axis, respectively, e.g., P_x>P_y and S_x>S_y; otherwise, such as P_x>P_y and , the negative anisotropy.

2.2 Speed and acceleration dynamics

The migration trajectories can be generated by APRW model with certain motile parameters which are set by referencing experimental data [33, 62]. In all simulations in the paper, the time step \text{d}t=0.01 min, and the recording time interval \text{Δ}t=1 min. Each trajectory is recorded with a length of 10⁴ min, giving the total number of samples N=10⁴.

Figure 1(a) gives a single trajectory generated by APRW model. With the trajectory r(t), the velocity vector v(t) can be calculated with

\begin{array}{}\text{(10)}& v(t)=\frac{r(t)-r(t-\text{Δ}t)}{\text{Δ}t}\mathrm{.}\end{array}

The corresponding speeds v(t)=\sqrt{v_x^2+v_y^2} of individual migration trajectory are plotted in Fig. 1(b). Figure 1(c) plots 1000 trajectories exhibiting clearly an anisotropic distribution. With a larger memory time of P_x=20 min, the cell migration typically has less stochasticity along the x axis, so the distance explored by cell along the x axis is larger than that along the y axis. The averaged speed \langle v(t)\rangle over 1000 trajectories is shown as a function of time in Fig. 1(d).

For PRW model described by Eq. (1), the expected value of the acceleration a_v for a given speed v can be given as [30]

\begin{array}{}\text{(11)}& {\langle a_v\rangle }_v={\langle \frac{\text{d}v}{\text{d}t}\rangle }_v=-\frac{1}{P}v\mathrm{.}\end{array}

For APRW model, the similar relationships are held for acceleration components at x and y axes, respectively. As a result, we have

\begin{array}{}\text{(12)}& {\langle a_v\rangle }_v={\langle \frac{\text{d}v}{\text{d}t}\rangle }_v=-\frac{1}{P_x}{v}_x-\frac{1}{P_y}{v}_y\mathrm{.}\end{array}

In simulation, for any two successive velocity vectors of each migration trajectory, we first compute the components of the second velocity vector which are parallel and orthogonal to the first velocity vector, and then calculate the acceleration components which are parallel and orthogonal (i.e., a_p and a_o) to the first velocity vector, respectively. Based on the relationship of parallel acceleration a_p(t) and speed v(t) against time, one can obtain the average parallel acceleration \langle a_p\rangle for each given speed v, giving a speed-acceleration profile.

Figure 2 discusses the influence of the heterogeneity of persistence time P and migration speed S on the profiles between speed v and parallel acceleration \langle a_p\rangle. As expected, for the case of no heterogeneity of P and S, i.e., {\text{Φ}}_P=1 and {\text{Φ}}_S=1, the parallel acceleration decreases linearly with increasing speed. With the increase of {\text{Φ}}_P but still keeping {\text{Φ}}_S=1, the v–\langle a_p\rangle profiles become more and more nonlinear. However, as shown by Eq. (12), the change in {\text{Φ}}_S has little effect on the v–\langle a_p\rangle profiles when {\text{Φ}}_P=1, as plotted in Fig. 2(b).

Figures 3(a) and (b) show the influence of heterogeneity of persistence time (P_x\mathrm{,}{P}_y) on the v–\langle a_p\rangle profiles in the cases of positive and negative anisotropy, respectively. With the increase of persistence time P_y, the v–\langle a_p\rangle profiles gradually become nonlinear, giving a concave curve for positive anisotropy and a convex curve for negative anisotropy, respectively. As a result, according to the shape of speed-acceleration profile, one can distinguish the positive anisotropy from the negative anisotropy.

Figures 3(c) and (d) show the influence of the heterogeneity of migration speed on the v–\langle a_p\rangle profiles for positive and negative anisotropy, respectively. With the decrease of S_y or S_x(3.0–0.5 μm/min), the v–\langle a_p\rangle profiles move to lower right corner with a larger abscissa and a smaller ordinate, respectively.

2.3 Cell migration kinetics

2.3.1 Singular vector decomposition

As a fact, the migration directions of P_x and P_y of cells in experiment are generally unknown. Thus, as suggested in Refs. [57-59], we apply the singular vector decomposition (SVD) to define the intrinsic axes of the migration trajectory. With SVD, the velocity matrix M for a trajectory can be expressed as

\begin{array}{}\text{(13)}& M=U\Sigma V^{\ast }\mathrm{,}\end{array}

where U is the matrix eigenvector of the product M{M}^{\ast }, V is the matrix eigenvectors of the product M^{\ast }M, \Sigma is the singular values of the matrix M, and * denotes the transposed matrix. The first and second eigenvectors of V^{\ast } indicate the primary direction (p axis) and non-primary direction (np axis) of the trajectory, respectively. For a given migration trajectory comprising N coordinate points in 2D microenvironment, the matrix M is N\times 2 dimensional with the first column of the matrix containing the abscissa values of N points and the second column containing the ordinate values. The computed matrix V^{\ast } is 2\times 2 dimensional in which the first column represents the primary vector p and the second column represents the non-primary vector np. Thus, one can define the intrinsic orthogonal coordinates of p and np axes for each trajectory.

2.3.2 The velocity power spectrum

With the intrinsic orthogonal coordinates processed for each trajectory individually, the velocity components on primary and non-primary axes can be calculated as a function of time. With the Fourier-transform, the speed power spectra on primary and non-primary axes can be obtained for each trajectory individually. Finally, the averaged speed power spectra on primary and non-primary axes can be calculated with 1000 trajectories for a given APRW model. In the following, we propose that the four motile parameters of APRW model can be estimated by fitting the averaged speed power spectrum to the theoretical speed power spectrum.

It has been suggested that the power spectrum of velocity of cell migration with OU process can be given by [48]

\begin{array}{}\text{(14)}& P_u(f_k)=P_u^{(true)}(f_k)+\frac{4{\sigma }_{pos}^2}{\text{Δ}t}[1-\cos (\pi f_k/f_{Nyq})]\mathrm{,}\end{array}

where the first term on the right side of Eq. (14) is the true expression of power spectrum, with the following form

\begin{array}{}\text{(15)}& \begin{array}{ll}{P}_u^{(true)}(f_k)=\hfill & \frac{(1-c^2)}{c}{\left(\frac{P}{\text{Δ}t}\right)}^2{P}_v^{(aliased)}(f_k)\hfill \\ \hfill & +4D\left(1-\frac{1-c^2}{2c}\frac{P}{\text{Δ}t}\right)\mathrm{,}\hfill \end{array}\end{array}

in which the term P_v^{(aliased)}(f_k) is defined by

\begin{array}{}\text{(16)}& P_v^{(aliased)}(f_k)\equiv \frac{\langle |\widehat{v}{|}^2\rangle }{t_{msr}}=\frac{(1-c^2)2D\text{Δ}t/P}{1+c^2-2c\cos (\pi f_k/f_{Nyq})}\mathrm{,}\end{array}

and the second term on the right side of Eq. (14) is an additional noise term when considering the effect of localization noise. Here, c=\exp (-\text{Δ}t/P), f_k=k\cdot \text{Δ}f, \text{Δ}f=1/t_{msr}, t_{msr}=N\cdot \text{Δ}t and f_{Nyq}=1/(2\text{Δ}t). N is the total number of recording time, \text{Δ}t is the sampling time interval with N=10⁴, and \text{Δ}t=1 min here. Considering the absence of localization noise in computer simulation, the noise parameter {\sigma }_{pos} is set as 0.0001 in our simulation.

As a fact, the trajectory of APRW process is the overlap of two independent PRW processes in x and y axes. Thus, the true expression of power spectrum of cell speed in APRW model should be the summation of two exponential decays, giving a double exponential decay, i.e.,

\begin{array}{}\text{(17)}& P_u(f_k)=P_u^1(f_k)+P_u^2(f_k)\mathrm{,}\end{array}

where P_u^1(f_k) and P_u^2(f_k) are given by Eq. (14) with two decay modes.

Such a double exponential decay for speed power spectrum will typically generate the double exponential decays for power spectra both in primary and non-primary directions if these directions are different from the original x and y axes of APRW model. As a result, our goal is to fit the theoretical expression of Eq. (17) with experimental power spectra in primary and non-primary directions, which will return a set of motile parameters \theta =(D_1\mathrm{,}{P}_1\mathrm{,}{D}_2\mathrm{,}{P}_2) in each direction.

In this paper, we apply the maximum likelihood estimation to obtain the motile parameters of \theta =(D_1\mathrm{,}{P}_1\mathrm{,}{D}_2\mathrm{,}{P}_2). For a given power spectrum {\left\{|\widehat{u}{|}^2/t_{msr}\right\}}_{k=\mathrm{0,}\cdots \mathrm{,}N-1}, the log-likelihood function is

\begin{array}{}\text{(18)}& \mathcal{l}\left(\theta |{\left\{\frac{|{\widehat{u}}_k{|}^2}{t_{msr}}\right\}}_{k=\mathrm{0,}\cdots \mathrm{,}N-1}\right)=2{\displaystyle \sum _{k=1}^N}\log \left(\frac{2}{P_u(f_k)}\right)+{\displaystyle \sum _{k=1}^N}\log \left(\frac{|{\widehat{u}}_k{|}^2}{t_{msr}}\right)-{\displaystyle \sum _{i=1}^N}\left(\frac{2}{P_u(f_k)}\cdot \frac{|{\widehat{u}}_k{|}^2}{t_{msr}}\right)\mathrm{,}\end{array}

where P_u(f_k) is theoretical power spectral value depending on APRW model with \theta =(D\mathrm{,}P\mathrm{,}{\sigma }_{pos}). The relationship between diffusion coefficient D and persistence time P is given as follows:

\begin{array}{}\text{(19)}& S=\sqrt{\frac{2D}{P}}\mathrm{.}\end{array}

As a result, the four motile parameters of P₁, P₂, S₁ and S₂ can be obtained with such a fitting method for each power spectrum in primary or non-primary direction.

3 Results

Now we apply our method to derive the motile parameters of APRW trajectories. In the simulation, we consider different migration models with different combinations of motile parameters, including high heterogeneity and low heterogeneity, as shown in Table 1.

The power spectra for APRW model with high and low heterogeneities are shown in Figs. 4(a) and (b), respectively.

As shown in Fig. 4, the power spectra can be divided into two parts: the horizontal part at low frequency range and attenuating part at high frequency range. The horizontal part denotes the power spectrum of Gaussian white noise, while the attenuating part manifests the persistence migration of cells which can be characterized by P and S with Eq. (17).

For a large HI in Fig. 4(a), the values of power spectrum of the horizontal part in p axis are larger than those in np axis, denoting that the value of diffusion coefficient D in p axis is larger than that in np axis. The nonlinearity of the attenuating part reveals the anisotropy of cell migration of APRW model. For a small HI in Fig. 4(b), the overlapping of the power spectra indicates that the related motile parameters may be similar in both directions.

With Eq. (17), the attenuating part of power spectrum is fitted with two exponential decays in each axis. As a result, we have two sets of fitted parameters of P and S in each direction. Figure 5(a) plots the two larger fitted P in p and np directions, and 5(b) plots two smaller P, as a function of theoretical P. Then, the final two persistence times are defined as the average of the two large and two small fitted persistence times, respectively, which are shown with triangles in Figs. 5(a) and (b), respectively.

Figures 5(c) and (d) plot the corresponding migration speeds S calculated by Eq. (19), respectively, as a function of theoretical value of S. Obviously, the fitted results of migration speeds all give certain deviations from the theoretical values. As a result, we furthermore consider an adjustment for the parameter of migration speeds, which is given as

\begin{array}{}\text{(20)}& S_p=\sqrt{S_{p1}^2+{[{\gamma }_p\cdot (S_{p2}+S_{np1})]}^2}\mathrm{,}\end{array}

\begin{array}{}\text{(21)}& S_{np}=\sqrt{S_{np2}^2+{[{\gamma }_{np}\cdot (S_{p2}+S_{np1})]}^2}\mathrm{.}\end{array}

The adjusting formulas are for migration speeds in p and np axes, respectively. Here, {\gamma }_p and {\gamma }_{np} are two adjusting factors. For positive anisotropy in Figs. 5(c) and (d), {\gamma }_p=\max \{\gamma \mathrm{,1}-\gamma \} and {\gamma }_{np}=\min \{\gamma \mathrm{,1}-\gamma \}, and for negative anisotropy in Figs. 6(c) and (d), {\gamma }_p=\min \{\gamma \mathrm{,1}-\gamma \} and {\gamma }_{np}=\max \{\gamma \mathrm{,1}-\gamma \}, in which γ is a weighing factor given by

\begin{array}{}\text{(22)}& \gamma =\frac{S_{np1}^2}{S_{np1}^2+S_{np2}^2}\mathrm{.}\end{array}

Note that the value of γ in each direction differs from each other. The corrected migration speeds are shown in green triangle in Figs. 5(c) and (d), presenting a better fitting to the theoretical values.

To discuss the universality of the method in dealing with different migration heterogeneity, we also analyzed the case with a low heterogeneity. The calculation results are shown in Fig. 6. Figures 6(a) and (b) show that there are certain deviations between the fitted persistence times and theoretical values both in p and np axes. However, the average of two fitted persistence times is closer to the theoretical value, which is shown as the triangles in Figs. 6(a) and (b). The calculated migration speeds are plotted in Figs. 6(c) and (d), in which the corrected S with Eqs. (20) and (21) are well close to the theoretical values indicating the generality and effectiveness of the adjusting method for migration speed.

4 Conclusions and discussion

Cell migration is a critical process in the development and maintenance of multicellular organisms. The cell migration in anisotropic microenvironment can be described by APRW model which is characterized by four motile parameters, i.e., two persistence times and two migration speeds at perpendicular directions. Thus, it is of importance to derive these motile parameters from the migration trajectories, in order to characterize the features of cell migration. In this work, we develop a method to obtain motile parameters of APRW trajectories.

We first analyze the migration trajectory produced by APRW model and find a nonlinear relationship between speeds and accelerations of cell trajectories in the forward direction, which is different from the linear relationship corresponding to PRW model. Second, we show that one can distinguish the positive and negative anisotropy from the shape of speed-acceleration profile, which reflects different relationships between persistence time and migration speed. Then the singular vector decomposition is applied to each migration trajectory individually to define the intrinsic primary and non-primary axes. Furthermore, one can calculate the speed power spectra in the intrinsic primary and non-primary axes for each migration trajectory. We indicate that the motile parameters of APRW model can be fitted with the averaged speed power spectra to the theoretical expression based on the maximum likelihood estimation.

APRW process generates two different single exponential decays for speed power spectra in original x and y axes, respectively, giving typically the double exponential decays for speed power spectra in all other directions. Thus, if the calculated primary and non-primary directions are different from the original x and y axes of APRW model, which is usually the case, the averaged speed power spectra in calculated primary and non-primary directions should be fitted by double exponential decays.

As a result, fitting with the speed power spectrum either in primary direction or in non-primary direction, one can obtain the four motile parameters of APRW process. With two sets of fitted motile parameters, our simulation results indicate that the direct average of the two large persistence times and the two small persistence times can give reliable values of the true persistence times. However, the migration speeds have to be corrected according to Eqs. (20) and (21) in order to obtain a good estimation of the parameters.

In the paper, to characterize the moving behaviors of cells, we propose a method to derive the motile parameters of cell migration in anisotropic microenvironment. In our simulation, the fitting results are obtained from an average of a large number of computer trajectories, which may be different from the experimental situation where short trajectories are recorded for not so many cell trajectories. In detail, we consider a length of 10⁴ min with 10⁴ points for 1000 trajectories for a better estimation of motile parameters. In the typical experiment, trajectories with a length of about 300 min for less than 100 cells are recorded. Our simulations show that with such trajectories, one could not derive reliable motile parameters. Furthermore, various noises are always presented in experiment, which is not considered in the simulation. Thus, how to improve the method to deal with the biological experimental data of cell migration is the future research interest.

CHANGE HISTORY

A Correction to this paper has been published:
https://doi.org

The Retraction Note was updated to reflect reception of statements of disagreement from all authors.

This article has been retracted.Please see the Retraction Notice for more detail:
https://doi.org

References

[1] M. Vicente-Manzanares and A. R. Horwitz, Cell migration: An overview, Methods Mol. Biol. 769, 1 (2011)

[2] D. S. Vasilev, N. M. Dubrovskaya, N. L. Tumanova, and I. A. Zhuravin, Prenatal hypoxia in different periods of embryogenesis differentially affects cell migration, neuronal plasticity, and rat behavior in postnatal ontogenesis, Front. Neurosci. 10, 126 (2016)

[3] M. Zalokar and I. Erk, J. Microsc. Biol. Cell. 25, 97 (1976)

[4] P. Kulesa, D. L. Ellies, and P. A. Trainor, Comparative analysis of neural crest cell death, migration, and function during vertebrate embryogenesis, Dev. Dyn. 229(1), 14 (2004)

[5] W. S. Krawczyk, A pattern of epidermal cell migration during wound healing, J. Cell Biol. 49(2), 247 (1971)

[6] G. D. Sharma, J. He, and H. E. P. Bazan, p38 and ERK1/2 coordinate cellular migration and proliferation in epithelial wound healing: Evidence of cross-talk activation between MAP kinase cascades, J. Biol. Chem. 278(24), 21989 (2003)

[7] T. Tao, A. Robichaud, S. Nadeau, R. Savoie, B. Gallant, and R. Ouellette, Effect of cumulus cell removal on the fertilization and the day 3 embryo quality in human IVF, International Congress Series 1271, 135 (2004)

[8] J. Bowszyc, J. Bowszyc, and T. Machońko, 15-years of studies of the immunological cellular response in syphilis in humans using migration inhibition tests, Przegl. Dermatol. 72(2), 134 (1985)

[9] H. Zhang, Y. Han, J. Tao, S. Liu, C. Yan, and S. Li, Cellular repressor of E1A-stimulated genes regulates vascular endothelial cell migration by the ILK/AKT/mTOR/VEGF(165) signaling pathway, Exp. Cell Res. 317(20), 2904 (2011)

[10] X. L. Fang Wei, M. Cai, Y. Liu, P. Jung, and J. W. Shuai, Regulation of 1, 4, 5-triphosphate receptor channel gating dynamics by mutant presenilin in Alzheimer’s disease cells, Front. Phys. 12(3), 128702 (2017)

[11] S. X. Liu, Y. Z. Geng, and S. W. Yan, Structural effects and competition mechanisms targeting the interactions between p53 and MDM2 for cancer therapy, Front. Phys. 12(3), 128908 (2017)

[12] H. C. Berg and F. Dyson, Random walks in biology, Phys. Today 40(3), 73 (1987)

[13] L. X. Li, Photon diffusion in a relativistically expanding sphere, Front. Phys. 8(5), 555 (2013)

[14] X. H. Li, G. Yang, and J. P. Huang, Chaotic-periodic transition in a two-sided minority game, Front. Phys. 11(4), 118901 (2016)

[15] S. Huang, C. P. Brangwynne, K. K. Parker, and D. E. Ingber, Symmetry-breaking in mammalian cell cohort migration during tissue pattern formation: Role of random-walk persistence, Cell Motil. Cytoskeleton 61(4), 201 (2005)

[16] Z. X. Niu, T. Hang, and Y. Chen, Molecular dynamics study of nanodroplet diffusion on smooth solid surfaces, Front. Phys. 13(5), 137804 (2018)

[17] Y. A. Yan and J. S. Shao, Stochastic description of quantum Brownian dynamics, Front. Phys. 11(4), 110309 (2016)

[18] M. Boguñá, J. M. Porrà, and J. Masoliver, Generalization of the persistent random walk to dimensions greater than 1, Phys. Rev. E 58(6), 6992 (1998)

[19] C. Peggy, A. L. Ny, B. D. Loynes, et al., Persistent random walks (I): Recurrence versus transience, J. Theor. Probab. 31(1), 1 (2016)

[20] H. I. Wu, B. L. Li, T. A. Springer, and W. H. Neill, Modelling animal movement as a persistent random walk in two dimensions: Expected magnitude of net displacement, Ecol. Modell. 132(1–2), 115 (2000)

[21] M. Schienbein and H. Gruler, Langevin equation, Fokker-Planck equation and cell migration, Bull. Math. Biol. 55(3), 585 (1993)

[22] D. S. Lemons and A. Gythiel, Paul Langevin’s 1908 paper “On the Theory of Brownian Motion” [Sur la théorie du mouvement brownien, C. R. Acad. Sci. (Paris) 146, 530–533 (1908)], Am. J. Phys. 65(11), 1079 (1997)

[23] J. L. Parada, G. Carrillo-Castañeda, and M. V. Ortega, Profile of the enzymes of the Krebs cycle in Salmonella typhimurium during the utilization of succinate, acetate, and citrate for growth, Rev. Latinoam. Microbiol. 15, 29 (1973)

[24] C. S. Patlak, Random walk with persistence and external bias, Bull. Math. Biophys. 15(3), 311 (1953)

[25] B. Wang, J. Kuo, S. C. Bae, and S. Granick, When Brownian diffusion is not Gaussian, Nat. Mater. 11(6), 481 (2012)

[26] T. H. Harris, E. J. Banigan, D. A. Christian, C. Konradt, E. D. Tait Wojno, K. Norose, E. H. Wilson, B. John, W. Weninger, A. D. Luster, A. J. Liu, and C. A. Hunter, Generalized Lévy walks and the role of chemokines in migration of effector CD8+ T cells, Nature 486(7404), 545 (2012)

[27] D. W. Sims, E. J. Southall, N. E. Humphries, G. C. Hays, C. J. A. Bradshaw, J. W. Pitchford, A. James, M. Z. Ahmed, A. S. Brierley, M. A. Hindell, D. Morritt, M. K. Musyl, D. Righton, E. L. C. Shepard, V. J. Wearmouth, R. P. Wilson, M. J. Witt, and J. D. Metcalfe, Scaling laws of marine predator search behaviour, Nature 451(7182), 1098 (2008)

[28] A. Czirók, K. Schlett, E. Madarasz, and T. Vicsek, Exponential distribution of locomotion activity in cell cultures, Phys. Rev. Lett. 81(14), 3038 (1998)

[29] A. Upadhyaya, J. P. Rieu, J. A. Glazier, and Y. Sawada, Anomalous diffusion and non-Gaussian velocity distribution of Hydra cells in cellular aggregates, Physica A 293(3–4), 549 (2001)

[30] D. Selmeczi, S. Mosler, P. H. Hagedorn, N. B. Larsen, and H. Flyvbjerg, Cell motility as persistent random motion: Theories from experiments, Biophys. J. 89(2), 912 (2005)

[31] L. Liu, G. Duclos, B. Sun, J. Lee, A. Wu, Y. Kam, E. D. Sontag, H. A. Stone, J. C. Sturm, R. A. Gatenby, and R. H. Austin, Minimization of thermodynamic costs in cancer cell invasion, Proc. Natl. Acad. Sci. USA 110(5), 1686 (2013)

[32] J. Zhu, L. Liang, Y. Jiao, and L. Liu, Enhanced invasion of metastatic cancer cells via extracellular matrix interface, PLOS One 10(2), e0118058 (2015)

[33] P. H. Wu, A. Giri, S. X. Sun, and D. Wirtz, Three-dimensional cell migration does not follow a random walk, Proc. Natl. Acad. Sci. USA 111(11), 3949 (2014)

[34] P. H. Wu, A. Giri, and D. Wirtz, Statistical analysis of cell migration in 3D using the anisotropic persistent random walk model, Nat. Protoc. 10(3), 517 (2015)

[35] B. B. Mandelbrot and J. W. Van Ness, Fractional Brownian motions, fractional noises and applications, SIAM Rev. 10(4), 422 (1968)

[36] R. Metzler and J. Klafter, The random walk’s guide to anomalous diffusion: A fractional dynamics approach, Phys. Rep. 339(1), 1 (2000)

[37] B. Liu and J. Goree, Superdiffusion and non-Gaussian statistics in a driven-dissipative 2D dusty plasma, Phys. Rev. Lett. 100(5), 055003 (2008)

[38] M. Weiss, M. Elsner, F. Kartberg, and T. Nilsson, Anomalous subdiffusion is a measure for cytoplasmic crowding in living cells, Biophys. J. 87(5), 3518 (2004)

[39] C. L. Stokes, D. A. Lauffenburger, and S. K. Williams, Migration of individual microvessel endothelial cells: Stochastic model and parameter measurement, J. Cell. Sci. 99(Pt 2), 419 (1991)

[40] L. Li, E. C. Cox, and H. Flyvbjerg, “Dicty dynamics”: Dictyostelium motility as persistent random motion, Phys. Biol. 8(4), 046006 (2011)

[41] F. Höfling and T. Franosch, Anomalous transport in the crowded world of biological cells, Rep. Prog. Phys. 76(4), 046602 (2013)

[42] V. Zaburdaev, S. Denisov, and J. Klafter, Lévy walks, Rev. Mod. Phys. 87(2), 483 (2015)

[43] R. T. Tranquillo, D. A. Lauffenburger, and S. H. Zigmond, A stochastic model for leukocyte random motility and chemotaxis based on receptor binding fluctuations, J. Cell Biol. 106(2), 303 (1988)

[44] G. A. Dunn, Characterising a kinesis response: Time averaged measures of cell speed and directional persistence, Agents Actions Suppl. 12, 14 (1983)

[45] C. L. Stokes and D. A. Lauffenburger, Analysis of the roles of microvessel endothelial cell random motility and chemotaxis in angiogenesis, J. Theor. Biol. 152(3), 377 (1991)

[46] M. R. Parkhurst and W. M. Saltzman, Quantification of human neutrophil motility in three-dimensional collagen gels: Effect of collagen concentration, Biophys. J. 61(2), 306 (1992)

[47] R. Gorelik and A. Gautreau, Quantitative and unbiased analysis of directional persistence in cell migration, Nat. Protoc. 9(8), 1931 (2014)

[48] J. N. Pedersen, L. Li, C. Grǎdinaru, R. H. Austin, E. C. Cox, and H. Flyvbjerg, How to connect time-lapse recorded trajectories of motile microorganisms with dynamical models in continuous time, Phys. Rev. E 94(6–1), 062401 (2016)

[49] A. P. Dempster, N. M. Laird, and D. B. Rubin, Maximum likelihood from incomplete data via the EM algorithm, J. R. Stat. Soc. B 39(1), 1 (1977)

[50] C. L. Vestergaard, J. N. Pedersen, K. I. Mortensen, and H. Flyvbjerg, Estimation of motility parameters from trajectory data, Eur. Phys. J. Spec. Top. 224(7), 1151 (2015)

[51] C. L. Vestergaard, P. C. Blainey, and H. Flyvbjerg, Optimal estimation of diffusion coefficients from single-particle trajectories, Phys. Rev. E 89(2), 022726 (2014)

[52] J. H. Chen and H. Y. Fan, On the core of the fractional Fourier transform and its role in composing complex fractional Fourier transformations and Fresnel transformations, Front. Phys. 10(1), 100301 (2015)

[53] Z. K. Wu, P. Li, and Y. Z. Gu, Propagation dynamics of finite-energy Airy beams in nonlocal nonlinear media, Front. Phys. 12(5), 124203 (2017)

[54] K. Katoh, K. Misawa, K. Kuma, and T. Miyata, MAFFT: A novel method for rapid multiple sequence alignment based on fast Fourier transform, Nucleic Acids Res. 30(14), 3059 (2002)

[55] S. B. Weinstein and P. M. Ebert, Data transmission by frequency division multiplexing using the discrete Fourier transform, IEEE Trans. Commun. Tech. 19, 628 (1971)

[56] S. F. Nørrelykke and H. Flyvbjerg, Power spectrum analysis with least-squares fitting: Amplitude bias and its elimination, with application to optical tweezers and atomic force microscope cantilevers, Rev. Sci. Instrum. 81(7), 075103 (2010)

[57] Y. Liu, X. Zhang, Y. Wu, W. Liu, X. Li, R. Liu, L. Liu, and J. Shuai, Derivation of persistent time for anisotropic migration of cells, Chin. Phys. B 26(12), 128707 (2017)

[58] N. Le Bihan and J. Mars, Singular value decomposition of quaternion matrices: A new tool for vector-sensor signal processing, Signal Processing 84(7), 1177 (2004)

[59] M. E. Wall, A. Rechtsteiner, and L. M. Rocha, Singular value decomposition and principal component analysis, arXiv: physics/0208101v4 (2002)

[60] K. Berg-Sørensen and H. Flyvbjerg, Power spectrum analysis for optical tweezers, Rev. Sci. Instrum. 75(3), 594 (2004)

[61] J. R. Xie and B. H. Wang, Modularity-like objective function in annotated networks, Front. Phys. 12(6), 128903 (2017)

[62] R. Gorelik and A. Gautreau, The Arp2/3 inhibitory protein arpin induces cell turning by pausing cell migration, Cytoskeleton (Hoboken) 72(7), 362 (2015)

Rights & permissions

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature 2020