Skip to main content
Advertisement
  • Loading metrics

Anomalous diffusion and asymmetric tempering memory in neutrophil chemotaxis

  • Peter Dieterich ,

    Roles Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing

    peter.dieterich@tu-dresden.de

    Affiliation Institut für Physiologie, TU Dresden, Dresden, Germany

  • Otto Lindemann,

    Roles Data curation, Investigation, Methodology, Resources, Writing – review & editing

    Affiliation Institut für Physiologie II, Westfälische Wilhelms-Universität Münster, Münster, Germany

  • Mats Leif Moskopp,

    Roles Data curation, Formal analysis, Investigation, Methodology, Resources, Software, Writing – review & editing

    Affiliations Institut für Physiologie, TU Dresden, Dresden, Germany, Klinik für Neurochirurgie, Vivantes Klinikum im Friedrichshain, Berlin, Germany

  • Sebastien Tauzin,

    Roles Conceptualization, Data curation, Investigation, Methodology, Resources, Writing – review & editing

    Affiliation Department of Biology, Utah Valley University, Orem, Utah, United States of America

  • Anna Huttenlocher,

    Roles Conceptualization, Investigation, Methodology, Resources, Writing – review & editing

    Affiliation Huttenlocher Lab, Department of Medical Microbiology, University of Wisconsin-Madison, Madison, Wisconsin, United States of America

  • Rainer Klages,

    Roles Conceptualization, Funding acquisition, Investigation, Methodology, Validation, Writing – review & editing

    Affiliations School of Mathematical Sciences, Queen Mary University of London, London, United Kingdom, Max Planck Institut für Physik komplexer Systeme, Dresden, Germany

  • Aleksei Chechkin,

    Roles Conceptualization, Funding acquisition, Investigation, Methodology, Validation, Writing – review & editing

    Affiliations Institute of Physics and Astronomy, University of Potsdam, Potsdam-Golm, Germany, Faculty of Pure and Applied Mathematics, Hugo Steinhaus Center, Wrocław University of Science and Technology, Wrocław, Poland, Institute for Theoretical Physics, NSC KIPT, Kharkov, Ukraine

  • Albrecht Schwab

    Roles Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing

    Affiliation Institut für Physiologie II, Westfälische Wilhelms-Universität Münster, Münster, Germany

Abstract

The motility of neutrophils and their ability to sense and to react to chemoattractants in their environment are of central importance for the innate immunity. Neutrophils are guided towards sites of inflammation following the activation of G-protein coupled chemoattractant receptors such as CXCR2 whose signaling strongly depends on the activity of Ca2+ permeable TRPC6 channels. It is the aim of this study to analyze data sets obtained in vitro (murine neutrophils) and in vivo (zebrafish neutrophils) with a stochastic mathematical model to gain deeper insight into the underlying mechanisms. The model is based on the analysis of trajectories of individual neutrophils. Bayesian data analysis, including the covariances of positions for fractional Brownian motion as well as for exponentially and power-law tempered model variants, allows the estimation of parameters and model selection. Our model-based analysis reveals that wildtype neutrophils show pure superdiffusive fractional Brownian motion. This so-called anomalous dynamics is characterized by temporal long-range correlations for the movement into the direction of the chemotactic CXCL1 gradient. Pure superdiffusion is absent vertically to this gradient. This points to an asymmetric ‘memory’ of the migratory machinery, which is found both in vitro and in vivo. CXCR2 blockade and TRPC6-knockout cause tempering of temporal correlations in the chemotactic gradient. This can be interpreted as a progressive loss of memory, which leads to a marked reduction of chemotaxis and search efficiency of neutrophils. In summary, our findings indicate that spatially differential regulation of anomalous dynamics appears to play a central role in guiding efficient chemotactic behavior.

Author summary

Neutrophil granulocytes are essential for the first host defense. After leaving the blood circulation they migrate efficiently towards sites of inflammation. They are guided by chemoattractants released from cells within the inflammatory foci. On a cellular level, directional migration is a consequence of cellular front-rear asymmetry which is induced by the concentration gradient of the chemoattractants. The generation and maintenance of this asymmetry, however, is not yet fully understood. Here we analyzed the paths of chemotacting neutrophils with different stochastic models to gain further insight into the underlying mechanisms. Wildtype chemotacting neutrophils show an anomalous superdiffusive behavior. CXCR2 blockade and TRPC6-knockout cause the tempering of temporal correlations and a reduction of chemotaxis. Importantly, such tempering is found both in vitro and in vivo. These findings indicate that the maintenance of anomalous dynamics is crucial for chemotactic behavior and the search efficiency of neutrophils.

Introduction

Being part of the innate immune system neutrophils play a crucial role during the first host defense. Their recruitment to a site of inflammation is an early event in the immune response and follows a well defined sequence of dynamical processes: rolling, adhesion, strengthening of adhesion, intravascular crawling and endothelial transmigration [1, 2]. After leaving the blood vessels neutrophils are guided by chemoattractant gradients to reach the site of inflammation [3]. On a molecular level chemotactic sensing involves the activation of G-protein coupled chemoattractant receptors such as CXCR2 whose activation trigger cytosolic signaling cascades underlying cell polarization in the chemotactic gradient and directed migration [4]. Chemotactic signaling frequently leads to a rise of the intracellular Ca2+ concentration by Ca2+ release from cytosolic stores and Ca2+ influx via Ca2+ permeable ion channels in the plasma membrane [58] such as members of the family of transient receptor potential (TRP) channels ([911] and [12] for a review). Impairing Ca2+ signaling interferes with neutrophil recruitment and chemotaxis [11]. We showed that TRPC6 channels are required for CXCR2-mediated chemotaxis of murine neutrophils [11]. Based on a conventional analysis of chemotaxis experiments that comprises speed, translocation and chemotaxis index, TRPC6−/− neutrophils have a chemotaxis defect that appears quite similar to that caused by the inhibition of CXCR2 receptors in wildtype neutrophils. Thus, this standard analysis is not sensitive enough to distinguish between the different roles of CXCR2 receptors and TRPC6 channels or to propose a signaling hierarchy.

Different mathematical models have been suggested to explain the complex mechanisms underlying chemotaxis. At the level of a continuum description, examples include the polarized LEGI-BEN (local excitation global inhibition-biased excitable network) [13], excitable signal transduction networks [14], or pseudopod-centered models [15]. Chemotaxis was also analyzed within the context of Lévy walks for T cell migration in lymphoid tissue under in vivo conditions [16] or within the framework of Langevin equations for Dictyostelium [17]. Typical random walk models include temporal exponential correlations [1821] and possible heterogeneities of cells [22]. Previously, we applied models of stochastic anomalous dynamics to analyze migration without chemotaxis and to quantitatively compare the phenotype of wildtype and mutant MDCK-F cells [23]. The latter are deficient for the sodium proton exchanger NHE1, a transport protein involved in persistent cell migration [2325]. We discovered that cells are migrating superdiffusively. Superdiffusion of cell migration can be caused by a temporal memory effect. This in turn, is due to the temporal correlations of successive steps of cell migration. Thus, superdiffusively migrating cells maintain their direction of movement for longer time periods as those performing an uncorrelated random walk. We also found that this superdiffusive mode of migration is not affected by NHE1 deficiency. NHE1 rather stabilizes the polarization of migrating cells. Such behavior of polarization has later been modeled as ‘polarized excitable network’ producing persistent migration in the absence of chemotactic stimulation [26]. Thus, the model-based mathematical analysis has provided a much deeper insight into the mechanisms of cell migration without influences of chemotactic stimuli.

We therefore set out to develop an extended mathematical model that is based on stochastic dynamics of individual cell paths in order to analyse chemotaxis of neutrophil granulocytes in more detail. We used the same data set as in our previous work [11] and complemented it with data from in vivo experiments performed in zebrafish embryos [27]. Our results show that chemotaxis (in vitro and in vivo) leads to an anomalous dynamics with reduced tempering in the direction of the chemotactic gradient. Tempering can be interpreted as progressive loss of memory of the migration machinery. In addition, we find a much stronger effect of CXCR2 blockade than TRPC6-knockout on anomalous behavior which translates into different search efficiencies of the neutrophils. Thus, the model-based analysis of the chemotactic behavior of neutrophils yields a depth of insight on the molecular and systems levels that could not have been achieved with conventional data analysis.

Results

Overview of chemotaxis experiments

Fig 1 summarizes the chemotaxis experiments with the trajectories of the murine neutrophil groups used in this study. Trajectories are normalized to common starting points with the coordinates (x, y) = (0, 0). Under control conditions wildtype neutrophils migrate efficiently along the CXCL1 (KC) gradient towards the positive y-axis (Fig 1A). This is indicated by the black line which depicts the mean position of the cell population during the course of the experiment. The black lines are much shorter in Fig 1B–1D because directed migration along the CXCL1 gradient is severely impaired by CXCR2 blockade (Fig 1B and 1D) and TRPC6 knock-out (Fig 1C and 1D).

thumbnail
Fig 1. Overview of trajectories and mean drift behavior for murine neutrophils.

(A)—(D) Paths of individual cells are normalized to a common starting point at the origin of the coordinate system. They are derived from different experimental conditions. Murine neutrophils were observed with a sampling interval of dt = 5 s over a time period of T = 30 min. Final cell positions are marked by colored circles (A-D). The black line indicates the mean position of the corresponding cell ensemble over the observation period T. The mean squared displacement of the cell ensemble after observation time T is represented as a red ellipse. The size of the ellipse indicates different propagation of cells in x- and y-direction. Murine neutrophils are chemotacting in a CLCX1 gradient along the positive y-axis. (A) Wildtype cells (WT). (B) Wildtype cells in the presence of the CXCR2 antagonist SB225002 (WT + SB). (C) TRPC6-knockout cells (TRPC6−/−). (D) TRCP6-knockout cells in the presence of CXCR2 antagonist SB225002 (TRPC6−/− + SB).

https://doi.org/10.1371/journal.pcbi.1010089.g001

The red ellipses represent the mean squared propagation of the cell positions in x- and y-direction with respect to their common center at the final time point of the experiments. The size of the ellipse is considerably smaller after CXCR2 blockade (Fig 1B). The asymmetry of the ellipse, which is particularly prominent for wildtype cells, indicates that murine neutrophils seem to migrate differently in x- and y-direction due to the influence of the chemotactic gradient. This asymmetry is lost when CRCX2 receptors are blocked or TRPC6 channels are absent.

In order to further analyze the chemotactic behavior we calculated the first and second moments of displacements from the starting positions for individual neutrophils as a function of time under different conditions according to Eqs 1619 (see Material and methods). The first moment characterizes directed motion (typically designated as drift), whereas the second moment contains additional information about the diffusive spreading of cells. The second moment is also called mean squared displacement. Fig 2A–2D show these moments for neutrophils from wildtype mice. Grey lines indicate the time-average results for each single cell (Eqs 16 and 17) and the green line shows the time- and ensemble-average of all cells as given by Eqs 18 and 19. As expected, there is a drift only in parallel to the chemotactic gradient (Fig 2C), which increases linearly in the direction of chemotaxis. It is absent perpendicular to the gradient (Fig 2A). Analogously we plotted the mean squared displacement in Fig 2B and 2D. Note that the second moment in the direction of chemotaxis is by more than one order of magnitude larger than in the perpendicular direction. This reflects the contribution of the drift. Fig 2E and 2F summarize the mean drift and mean squared displacement for all four experimental conditions studied in murine neutrophils. CXCR2 and TRPC6-knockout lead to an almost identical decrease of the mean drift velocity vd given as slope of 〈y(t)〉. Combined CXCR2 inhibition and TRCP6-knockout cause a small additional reduction. Wildtype cells have a drift velocity of ∼ 7.4 μm/min. Applying the CXCR2 inhibitor SB225002 reduces the drift to ∼ 25% of this control value. TRPC6-knockout cells also have a similarly diminished chemotatic drift, which is further reduced by blocking the CXCR2 receptor in these cells. TRPC6KO differentially affects drift and mean squared displacement: The drift is more strongly reduced by the TRPC6KO than the mean squared displacement. This is consistent with our previous observation [11] that TRPC6KO predominantly impairs the chemtotatic efficiency rather than the migration motor.

thumbnail
Fig 2. Overview of the first and second moments.

(A), (B) and (C), (D) show the first and second moments of wildtype neutrophils with respect to x- and y-direction, respectively. Thin grey lines indicate time-averaged moments of single cells according to Eqs 16 and 17. The thick green lines show the time-ensemble-average of the corresponding moment according to Eqs 18 and 19. (E) compares the mean first moment of different experimental conditions with murine neutrophils. (F) shows the differences of the second moments, i.e. the mean squared displacement for murine neutrophils.

https://doi.org/10.1371/journal.pcbi.1010089.g002

Velocity autocorrelations of chemotacting neutrophils

In order to analyze the migration behavior of the neutrophils in more detail we calculated the velocity autocorrelation functions based on Eqs 21 and 22. They provide information about the temporal memory of migrating neutrophils. Fig 3 shows the time-ensemble-averaged velocity autocorrelation functions 〈vcorr,x|y(t)〉 for murine neutrophils. In general, these velocity autocorrelations slowly decay as shown in Fig 3 indicating a certain kind of temporal persistence which cannot be described by an uncorrelated random walk process. The finding of strong temporal autocorrelations is consistent with the model of Andrew and Insall [28] who proposed that lamellipodia preferentially grow from already existing protrusions. The experimental velocity correlation functions show typical noisy scatterings for longer times (t ≳ 10 min) thereby impeding a direct quantification of the decay as purely exponential or power-law. The velocity autocorrelation contains the diffusive and drift components of the stochastic process: Starting from 〈v(t = 0)2〉 at time t = 0 it decays towards the mean squared drift velocity for long times. In addition, we analyzed the velocity cross-correlations 〈vcross,xy(t)〉 according to Eqs 23 and 24 and did not find couplings of velocities between x- and y-direction (see S3 Data).

thumbnail
Fig 3. Velocity autocorrelation functions.

The diagrams show the velocity autocorrelations in x- and y-direction. Thick colored lines indicate the corresponding time-ensemble-averaged velocity autocorrelation function 〈vcorr,x|y(t)〉 as defined in Eq 22 with repect to x- and y-direction in bright and dark colors, respectively. The results of the best theoretical model are shown as black lines together with a grey area of uncertainty (± one standard deviation). (A) Murine wildtype neutrophils. (B) Wildtype neutrophils with CXCR2 blocker SB225002. (C) TRPC6−/− neutrophils. (D) TRPC6−/− neutrophils with CXCR2 blocker SB225002.

https://doi.org/10.1371/journal.pcbi.1010089.g003

Development of an optimal stochastic model for neutrophil chemotaxis

The migratory behavior of neutrophils in a chemotactic gradient as shown in Fig 1 is governed by two main processes: The first process refers to the directed migration of the cell population reflecting the effect of the external force induced by the chemotactic gradient. This external force can be assumed as constant since the mean position of cells along the y-axis increases linearly in time (Fig 2E). Diffusive spreading constitutes the second process. However, the temporal persistence of trajectories questions whether the migration process follows an uncorrelated random walk. The slow temporal decay of the velocity autocorrelation function as shown in Fig 3 clearly indicates temporal correlations. This is beyond the characteristics of ordinary Brownian motion, which is a totally uncorrelated process [29, 30]. Thus, these observations suggest to model cell migration as a correlated stochastic process, where the position of the cell (x, y) is influenced by an external constant force modeled by vd (CXCL1 gradient) and an active stochastic force ζx|y(t) generated by the cellular migration apparatus: (1) This approach corresponds to an overdamped Langevin equation driven by (fractional) noise. Couplings of the motion in x- and y-direction can be neglected due to missing cross-correlations of velocities. As the experimental data of the velocity autocorrelation function do not allow a unique identification of its decay as exponential or power-law form, we tested the validity of three models (A—C) beyond the Ornstein-Uhlenbeck with exponential correlations (model D) [1822]. These models consider different types of temporal correlations and tempering mechanisms. Tempering accelerates the decay of autocorrelations. Thereby, it may allow a transition from anomalous to normal dynamics. The models are defined in the following.

Model A: Pure fractional Brownian noise with correlations [30, 31] (2) with the generalized diffusion coefficient DH and the Hurst coefficient H (0 < H < 1). Γ represents the gamma function. Model A is characterized by a temporal power-law decay ∼ τ2H−2 with τ = |tt′|. The model reduces to Brownian motion for H = 1/2 and shows superdiffusive behavior for H > 1/2. The special case H = 1/2 can be obtained from the asymptotic expansion of the gamma function for ϵ → 0 [32] and a representation of the delta function [33]. With ϵ = 2H − 1 the limit of ϵ → 0 corresponds to H → 1/2. Thus, the correlation function reduces to 〈v2(τ)〉A = 2 D1/2 δ(τ) for H = 1/2.

Model B: Fractional Brownian noise with exponential tempering [34] (3) with tempering time τ*. The addition of the exponential term to model B causes an exponential loss of correlations for times > τ*. Model B reduces to model A for τ* → ∞.

Model C: Fractional Brownian noise with power-law tempering [34] (4) with tempering time τ* and power-law tempering exponent μ ≥ 0. Model C is identical with model A for μ = 0 or τ* → ∞.

Model D: Ornstein-Uhlenbeck like processes are often used as reference processes to describe and quantify cell migration [1822]. Thus, we also included this approach, which is characterized by an exponentially decreasing velocity autocorrelation function: (5) Parameter denotes the mean squared velocity of cells and γ captures the decay of the exponential function with time-scale ∼ 1/γ.

Our models A—D are Gaussian processes that are fully characterized by their first and second moments. The first moments follow directly from Eq 1 and are identical for models A—D: (6) The second moments given by the covariances of these processes can be obtained by integrating the velocity autocorrelations of Eqs 25 twice: (7) They are different for each model: (8) (9) (10) (11) 1F1 and 2F1 represent the confluent hypergeometric and the Gaussian hypergeometric function, respectively [32]. The y-direction contains additional terms due to the drift velocity vd, where model parameters (H, DH, τ*, μ) and can take different values with respect to x- and y-direction. The corresponding second moments follow directly from the covariances of Eqs 811 for t = t′ as (12) (13) (14) (15) Analogously to the covariances, the msd in y-direction msdy(t)A|B|C|D is extended by the drift factor . In addition, the mean squared displacements for models A to C in Eqs 1214 reduce to normal diffusion 2 D1/2 t for H → 1/2 as 1F1(0, 2; z) = 1 and 2F1(μ, 0, 2; z) = 1 [32]. The msdx(t)D of the Ornstein-Uhlenbeck process in Eq 15 displays a ballistic phase for small times t ≪ 1/γ and normal diffusion for long times t ≫ 1/γ.

Application of the model to experimental data

Bayesian data analysis.

Bayesian data analysis [35, 36] is applied in order to link experimental data and the stochastic model for positions described in Eq 1. The Bayesian formalism including the covariances of the stochastic process (see Eqs 28 and 29) allows a logically consistent estimation of model parameters and their uncertainties. In addition, it provides the evidence for selecting the best model. Both calculations are based on the amount and quality of the experimental data. It is noteworthy, that the Bayesian formalism uses the positions of cell paths (raw experimental data) instead of performing a parameter fitting to the averaged mean squared displacement data. To avoid systematic errors due to pixelation of spatial coordinates, data were analyzed with dt = 10s, i.e. with every second data point.

Model classification.

In the first step we compared the four different models applied to all experimental data. Fig 4 summarizes the model probabilities for each experimental condition. We applied the models separately for x- and y-coordinates. It is clearly evident that movement of WT neutrophils along the chemotactic gradient (y-direction) is best modeled by pure fractional Brownian motion (model A). Whenever directional movement is absent (x-direction; Fig 4A) or impaired (y-direction of WT+SB, TRPC6KO and TRPC6KO+SB; Fig 4B), models with tempering are favored (model B and C). The necessity of introducing tempered models indicates that the extent of the temporal memory depends on the orientation with respect to the chemotactic gradient (x versus y) and on CXCR2 signaling and/or TRPC6 channels. Model probabilities of the Ornstein-Uhlenbeck process of model D were extremely low (< 6.510−51% for all experimental conditions) and therefore could not be made visible in Fig 4. Therefore, we omitted the Ornstein-Uhlenbeck process of model D from the further analysis. Nevertheless, the analysis revealed the time-scale 1/γ ∼ 0.1 … 0.2 min for all four experimental conditions. This indicates that a ballistic phase of migrating cells is hardly present and thus confirms the use of the overdamped modeling approach.

thumbnail
Fig 4. Model probabilities for murine neutrophils.

The model probabilities for the analyzed models A—C are color coded and shown with respect to (A) x- and (B) y-direction (probabilities of model D are too small to be visualized). Wildtype neutrophils follow pure fractional Brownian motion parallel to the chemotactic gradient, whereas the dynamics of neutrophils under all other conditions is either exponentially or power-law tempered.

https://doi.org/10.1371/journal.pcbi.1010089.g004

Quantitative analysis of model parameters.

To obtain an overview of the consequences of tempering we first analyzed the logarithmic derivative of the mean squared displacement βx|y(t) = d ln msdx|y(t)/d ln t (Fig 5A and 5B) perpendicular and parallel to the direction of the chemotactic gradient. Fig 5B illustrates that only WT neutrophils move with an untempered superdiffusive behavior with a constant βy ∼ 1.55 during the entire observation period. The parameter βy(t) of all other groups decays. The decay is more gradual in TRPC6KO neutrophils and stronger upon CXCR2 blockade. In contrast, the movement vertically to the chemotactic gradient (x-direction) is tempered for all four experimental groups in accordance with models B and C (Fig 5A).

thumbnail
Fig 5. Tempering behavior of murine neutrophils.

(A) and (B) show the temporal behavior of the logarithmic derivative β(t) of the mean squared displacement perpendicular (A) and parallel (B) to the direction of chemotaxis. The time-scale τ* of tempering is compared in (C) and (D) for the four different cell groups (τ* = ∞ for model A, pure fBm). (E) and (F) show the exponent μ of power-law tempering of model C. The color intensities of the parameter bars are proportional to the corresponding model probabilities given in Fig 4.

https://doi.org/10.1371/journal.pcbi.1010089.g005

A more quantitative analysis of the tempering behavior is obtained by considering the model parameters τ* (tempering time; model B/C) and μ (power-law tempering coefficient; model C). Numerical estimates and uncertainties of these parameters for all three models are plotted in Fig 5C–5F. The color intensities of the parameter bars are proportional to the corresponding model probabilities shown in Fig 4. The fact that model A (pure fractional Brownian motion) is favored for WT neutrophils is reflected by an infinite tempering time τ*. TRPC6KO neutrophils have a tempering time in the order of τ* ∼ 20 min. It is further decreased in the presence of the CXCR2 blocker.

The parameter μ describes the strength of power-law tempering in model C. Values of μ are between 1.5 and 3.0 except for WT neutrophils. Thus, the parameter μ is in the range of μ > 2H − 1 for WT+SB, TRPC6KO and TRPC6KO+SB groups. This indicates so-called strong power-law tempering [34]. Thus, TRPC6KO neutrophils and neutrophils with inhibited CXCR2 receptor perform normal diffusion for long times.

All three models contain the parameters H (Hurst coefficient), DH (diffusion constant) and vd (drift velocity) as shown in Fig 6. First, we consider the model parameters for movement in the y-direction in Fig 6B, 6D and 6F. Modeling recapitulates that drift velocity of neutrophils moving in a CXCL1 gradient strongly depends on CXCR2 signaling and the presence of TRPC6 channels. It is almost zero when the receptor is inhibited and / or the channel is absent. Model parameters vd are almost independent of the applied model A—C. The Hurst coefficient of WT neutrophils (Fig 6B) amounts to 0.76±0.01 for the preferred model A indicating a superdiffusive behavior over the entire observation interval. The Hurst coefficients of all other experimental groups are higher. However, superdiffusive dynamics for these cells is tempered as indicated by the preferred models B and C. Thus, superdiffusive behavior of WT+SB, TRPC6KO and TRPCKO+SB neutrophils is rapidly attenuated with time (see Fig 5). The generalized diffusion coefficient DH (Fig 6D) is about 39.01 ± 1.89 μm2/min for WT neutrophils. WT+SB and TRPC6KO neutrophils have a higher diffusion coefficient.

thumbnail
Fig 6. Comparison of model parameters for murine neutrophils.

Parameters are plotted in x- and y-direction together with their standard deviation for the three analyzed models: A—fractional Brownian motion (fBm), B—fBm with exponential tempering, and C—fBm with power-law tempering. The color intensities of the parameter bars are proportional to the corresponding model probabilities in Fig 4. (A) and (B) Hurst exponent H, (C) and (D) generalized diffusion coefficient, (E) and (F) absolute value of drift velocity vd.

https://doi.org/10.1371/journal.pcbi.1010089.g006

In summary, it is notable that parameters H and DH are quite similar in parallel and vertical to the chemotactic gradient (Fig 6). This could represent a common symmetric dynamics of the cells for small times. Our results also demonstrate that long term behavior is predominately determined by the degree of tempering (parameters τ* and μ). Asymmetric tempering is stronger vertical to the direction of chemotaxis. Less tempering in parallel to the chemotactic gradient may favor directed movement which is assessed by parameter vd.

Behavior of zebrafish neutrophils

So far we have analyzed the chemotactic behavior of murine neutrophils in vitro. In the next step we wanted to ensure that these results are also relevant for the in vivo situation. Therefore, we monitored directed movement of neutrophils in a zebrafish model. Neutrophil movement was induced by tail transection and monitored for up to 2 hours. Experiments delivered N = 108 and N = 106 trajectories with different lengths for control (WT) and SB225002-treated zebrafish (WT+SB), respectively. Fig 7 provides an overview of the trajectories obtained under both conditions (Fig 7A and 7C) and the resulting model parameters (Fig 7E–7J). It is evident that the model-based analysis of the trajectories of zebrafish neutrophils yields qualitatively almost the same results as described above for murine neutrophils (Figs 5 and 6).

thumbnail
Fig 7. Model parameters for the dynamics of neutrophils in zebrafish larvae.

(A) and (C) give an overview of all different trajectories of neutrophils moving in zebrafish larvae. Paths of neutrophils have different lengths. Neutrophils treated with blocker SB225002 (marked as WT+SB) have a reduced drift (C) compared to untreated cells (A). Model probabilities for WT neutrophils (B) and for the WT+SB group favor the exponentially tempered model under all directions and conditions (color coding as in Fig 4). (E) and (F) show the temporal behavior of the logarithmic derivative β(t) of the mean squared displacement for the two cell groups. WT+SB neutrophils are tempered much more strongly than WT neutrophils. The movement of WT neutrophils in parallel to the direction of chemotaxis (y-direction) is only weakly tempered. (G)—(J) show the values of model parameter for fBm with exponential tempering: (G) Hurst exponent, (H) diffusion coefficient DH, (I) absolute value of drift vd, and (J) tempering time τ*.

https://doi.org/10.1371/journal.pcbi.1010089.g007

Model probabilities for WT neutrophils (Fig 7B) and for the WT+SB group (Fig 7D) are clearly in favor of model B (fractional Brownian motion with exponential tempering). Only for WT neutrophils moving in the direction of chemotaxis, model A and model C show similar, but slightly smaller probabilities than model B. The effect of tempering is illustrated in Fig 7E and 7F by the decay of the logarithmic derivative β(t) of the mean squared displacement. Movements of WT neutrophils in y-direction are only weakly tempered so that their behavior almost fits to the pure fBm of model A. In contrast, tempering vertical to the chemotactic gradient (x-direction) is much stronger. As shown in Fig 7F, WT+SB neutrophils are tempered more strongly in both directions by the blockade of the CXCR2 receptor. Finally, parameters for neutrophils moving in zebrafish embryos show Hurst exponents H ∼ 0.80 − 0.86 (Fig 7G). The generalized diffusion coefficient DH is increased in parallel to the direction of chemotaxis for WT and WT+SB neutrophils (Fig 7H). This is most likely due to the anatomical structure of the tail fin of the zebrafish embryo [37]. It is characterized by rays that are aligned in parallel towards the the tip of the tail fin. This also applies to collagen fibers of the extracellular matrix and blood vessels [38]. Thus, these structural cues are in parallel to the direction of the chemotactic gradient which is generated by tail transection. They are relevant directional cues because the topography of the extracellular environment largely influences immune cell migration [39]. In addition to the increased diffusion coefficient, there is a clear directed drift velocity vd ∼ 2.5 μm/min, which is reduced to vd ∼ 0.8 μm/min by the blockade of the CXCR2 receptor (Fig 7I). Finally, tempering time τ* is around ∼ 25 min for WT cells in y-direction. This time is larger than than mean observation time of cells which is around ∼ 17.1 min and 19.5 min for WT and WT+SB, respectively.

Model-based simulation

Using the model parameters summarized in Figs 5 and 6 we simulated a large number of trajectories (100000) for each of our experimental conditions. The simulations were performed for the direction of chemotaxis (y-direction) using the corresponding covariance in Eqs 810 and parameters of the most probable model (Fig 4) for each experimental condition (see Material and methods for simulations based on the covariance function). I.e. model A with the corresponding WT-parameters were applied to simulate WT-conditions. Analogously, model C with the WT+SB parameters, model B with TRPC6KO parameters, and model C with TRPC6KO+SB parameters were used to simulate paths for WT+SB, TRPC6KO, and TRPC6KO+SB neutrophils, respectively.

We used these simulations to make predictions about the long-term search behavior of neutrophils and evaluated their search efficiency under the assumption that model parameters remain constant over the considered time period. The search behavior is quantified by the so-called first passage time tfpt describing the duration neutrophils need to pass a target for the first time [40]. Fig 8A shows the probability of first passage times P(tfpt) for neutrophils reaching a target at distance a = 200 μm apart from their starting point. This corresponds to the distance neutrophils are physiologically covering within the interstitium during swarming [41]. Wildtype neutrophils come in rapidly, whereas neutrophils from all other groups (TRPC6-knockout and/or CXCR2 blockade) arrive in a trickle. This is consistent with the kinetics of neutrophil swarming which begins after a lag phase of ∼ 15 min and which is impaired in CXCR2−/− neutrophils [42]. In Fig 8B the percentage of arriving cells after a time period of 5 h is plotted as a function of the target distance a. All simulated wildtype neutrophils reach the target within the observation period, while TRPC6-knock-out and CXCR2 blockade strongly reduce the number of arriving cells. Notably, the simulation results with parameters for TRPC6−/− neutrophils correspond quite well to our previous data on neutrophil recruitment in a peritonitis model [11]. Here, the number of recruited TRPC6−/− neutrophils is reduced to 60–70% of that in wildtype animals. The distance of capillaries to the peritonial surface is in the order of 100–400 μm [43]. Hence, such simulations could also be employed for planning the observation time interval e.g. when using end-point measurements of neutrophil emigration.

thumbnail
Fig 8. Comparison of arriving probabilities by simulations.

(A) Probability distributions of first passage times for simulated cells at a distance of 200 μm. Model parameters calculated for the respective experimental conditions were used for simulations as given in Figs 5 and 6. There is a surge of WT neutrophils reaching their target after 20 min, whereas all other cell types show a delayed arrival. (B) Percentage of cells arriving at the target that is located at a distance a from the starting points of the cell ensemble. All simulated wildtype cells arrive at the target within 5 h, whereas CXCR2 inhibition and/or TRPC6 ablation markedly reduce the number of cells arriving at their target.

https://doi.org/10.1371/journal.pcbi.1010089.g008

Discussion

We developed three stochastic models based on anomalous dynamics and tempering mechanisms and applied them together with the Ornstein-Uhlenbeck process to analyze neutrophil chemotaxis under different experimental conditions. Using a gradient of CXCL1 we examined chemotacting murine neutrophils under control conditions, in the presence of a CXCR2 blocker and in the absence of TRPC6 channels. In addition, murine neutrophils chemotacting in vitro were compared with the chemotactic behavior of neutrophils in zebrafish embryos. Our model-based approach reveals clear differences between the experimental conditions that were not as evident from previous conventional analysis [11]. CXCR2 inhibition has a more profound effect on the dynamics of chemotaxis than TRPC6-knockout suggesting that (i) TRPC6 channels are positioned downstream of the receptor and that (ii) CXCR2 receptors trigger other, TRPC6-independent intracellular signaling cascades. If all CXCR2 signaling had to pass via TRPC6 channels, its deletion would have elicited the identical phenotype. Yet, tempering and drift velocity are more dependent on CXCR2 receptors than on TRPC6 channels.

One of the key result is the finding that there is an asymmetry of anomalous migration dynamics in the direction of the chemotactic gradient and in transverse direction. Importantly, this asymmetric anomaly is also observed in vivo in the zebrafish embryo. This asymmetry appears to be a prerequisite for efficient chemotaxis. To the best of our knowledge this is the first description of asymmetric anomalous dynamics during chemotaxis in an intact organism as revealed by the strong and weak tempering behavior (see Fig 7E) in x- and y-direction, respectively. In an in vivo model of Toxoplasma gondii infection of the brain, migration of CD8+ T cells within the brain was described as a generalized Lévy walk, which was accelerated in a chemokinetic way by CXCL10 [16]. Thus, migration of CD8+ T cells appears to harbour anomalous properties as well [16]. Our analysis goes one step further in that we are analyzing chemotaxis and show the asymmetry and tempering of anomalous behavior parallel and perpendicular to the direction of the chemotactic gradient.

Moreover, we found that the anomalous properties of the neutrophils can be modulated by distinct molecular mechanisms. We discovered that the degree of anomaly and of its asymmetric tempering strongly depend on CXCR2 and to a lesser extent on TRPC6 channels. Accordingly, the parameter β(t) decays faster in blocker treated neutrophils than in TRPC6KO cells. As stated above, this also applies to neutrophils in zebrafish larvae. In our view it is likely that the higher degree of memory (corresponding to less tempering) in the direction of the chemotactic gradient favors the biased formation of pseudopods [15] or the biased activity of an excitable network and polarization underlying the polarized local excitation and global inhibition (LEGI-BEN) model [13]. Along these lines our analysis indicates that CXCR2 receptors and TRPC6 channels have a major impact on the bias of the excitable network as evidenced by the markedly reduced drift velocity when these membrane proteins are inhibited or absent. In contrast, TRPC6 channels have only a smaller effect on tempering and its asymmetry which could be reflected by the ‘polarization’ term in this model [26]. The combined effect, however, is sufficient to account for poor chemotaxis under both conditions. In addition, TRPC6 channels may be involved in the sensitivity of the gradient detection which depends on the polarization of the cells. In contrast, the CXCR2 inhibition disturbs drift velocity, tempering and asymmetry indicating that bias and polarization are strongly impaired.

Thus, by monitoring the migratory behavior on a ‘macroscopic’ level and performing a model-based analysis of the cell trajectories, we are able to obtain quantitative information on the molecular mechanisms underlying neutrophil chemotaxis. The effects of CXCR2 blockade or TRPC6-knockout may serve as proof-of-principle in this respect. As proposed earlier [44], CXCR2 receptors and TRPC6 channels seem to be signaling modules that transduce the ‘bias’ of newly formed pseudopodia [15] or of the excitable network within the lamellipodium [13]. Future studies will have to show whether this is linked to a polarized distribution of TRPC6 channel proteins or to localized activity of TRPC6 channels at the front of chemotaxing neutrophils. Then the channels could initiate local Ca2+ signaling domains as described earlier [6, 4548]. The latter possibility is supported by observations made in fibroblasts [49] where DAG [50], one of the main physiological activators of TRPC6 channels [51], was found to be enriched at the cell front.

The efficiency of the chemotaxis of neutrophils is related to the strength of tempering of the anomalous behavior: We observe optimal chemotaxis of WT neutrophils where tempering is absent in the direction of chemotaxis and present perpendicular to the chemotactic gradient. In this case anomality and temporal memory are maintained for the entire period of observation. In addition, the resulting drift velocity vd has a high value. Conversely, the absence of asymmetry and increased tempering in the direction of the gradient for TRPC6KO neutrophils and CXCR2 blockade is accompanied by an impaired chemotaxis, i.e. a reduced drift velocity vd. Then even the seemingly small additive effect of TRPC6 channel deletion and CXCR2 blockade [11] (see also Fig 2) translate to marked differences in search efficiency. For small times, the higher generalized diffusion coefficient DH of TRPC6KO versus WT+SB neutrophils favors the TRPC6KO first passage times. However, for longer times the higher drift velocity of WT+SB neutrophils gets dominant resulting in a higher percentage of arriving cells at a given target distance. These in silico results may have important implications for pathologies associated with CXCR2 signaling. One such example is the CXCR2-dependent recruitment of neutrophilic myeloid derived suppressor cells to cancers such as pancreatic ductal adenocarcinoma [52], breast cancer [53] or lung cancer [54]. In addition, our findings are also relevant for interpreting search behavior and optimal search strategies of anomalous processes [16, 5557].

Material and methods

Ethics statement

Experimental data of murine neutrophils are taken from Lindemann et al. [11] without further experiments (data from [11] were conducted with the permission 84–02.05.50.15.010 of the Landesamt für Natur, Umwelt und Verbraucherschutz Nordrhein-Westfalen, Germany). Adult and larval zebrafish were maintained according to protocols approved by the University of Wisconsin-Madison Institutional Animal Care and Use Committee [27].

Cells—Murine neutrophils

We further analyzed data sets of neutrophils [11] isolated from wildtype and TRPC6−/− mice migrating in chemotactic gradients of CXCL1 (keratinocyte-derived cytokine; CXCL1). Chemotaxis of both cell types was monitored in the absence and presence of the CXCR2 blocker SB225002 (500 nM). Chemotaxis experiments were performed in chemotaxis chambers (ibidi, Martinsried, Germany). Cells were seeded within a three-dimensional collagen I matrix (2.1 mg/ml; pH 7.4; BD, Heidelberg, Germany). Chemotaxis was recorded over a period of T = 30 min with a time sampling interval of dt = 5 s. Trajectories were obtained by tracking the movement of the cell center resulting in 361 x-y-positions for each cell path [11]. N = 60 and N = 100 tracks were used for the analysis of wild-type neutrophils and TRPC6−/− cells, respectively.

In vivo—Zebrafish neutrophils

In vivo experiments in previously described Tg (mpx:Dendra2) zebrafish larvae expressing green fluorescent neutrophils were performed to quantify neutrophil recruitment to wounds [27]. 3 dpf Tg (mpx:Dendra2) zebrafish larvae were pre-treated for 30 minutes with 5 μM of CXCR2 inhibitor SB225002 or with DMSO as control. Larvae were then tail transected and mounted in 1% agarose containing 5 μM SB225002 or DMSO. Neutrophil recruitment to the wound was recorded with a mean sampling interval of dt = 17 s up to 2 h with a laser scanning confocal microscope (FluoView FV1000, Olympus). The analysis of the video sequences and extraction of cell paths was performed with self-developed image processing software (based on python and the openCV library).

Data analysis

Quantitative data analysis and calculation of parameters were performed with extensions of programs developed by ourselves [23]. Cell trajectories were used to calculate the first and second moments. We analyze time-averages of single trajectories xk(t) of cell k as a function of time t in the interval [0, T]. These are defined as (16) (17) where the over-line indicates the time-average. The time-ensemble-average is given by the mean over all N cell paths k = 1 … N: (18) (19) The definitions for the y-direction apply analogously. The second moment is synonymously called mean squared displacement in the present work, i.e. msdx(t) = 〈x2(t)〉 and msdy(t) = 〈y2(t)〉. It should be noted here that the mean squared displacement is defined as non-central moment. We do not use central moments, because central moments in combination with time-ensemble averages and fractional dynamics cause several problems for systems with drift [58, 59]. In addition, our model analysis is based on the position of cells and not on the moments of positions. The velocity of migrating cells in x and y-direction vx|y(t) was calculated from trajectories as the difference quotient of two cell positions at times t + Δ and t. (20) For Δ the corresponding sampling interval of the experiments was applied. These velocities were used to calculate the autocorrelation function of velocities in x- and y-direction. Analogously to the time-average and time-ensemble-average definitions for the moments in Eqs 1619 velocity autocorrelations are defined as (21) (22) These definitions apply analogously to the y-direction. Accordingly, the time- and ensemble-averaged velocity cross-correlation function between velocities in x- and y-direction can be defined as (23) (24) to study possible couplings of velocities between different directions.

Bayesian inference of model parameters based on the covariance of the stochastic process

Gaussian stochastic processes as e.g. for positions x(t) are completely determined by their first moment μ(t) and covariance C: (25) An experimental sampling of a single trajectory of this process delivers M data points d = {d(t1), d(t2) … d(tM)} at times ti with i = 1 … M. The resulting likelihood for one cell path is given by [6063] (26) The mean drift μ and the covariance C may depend on model parameters θ with (27) (compare with Eqs 610). The likelihood of a cell ensemble consisting of N paths with data dk (k = 1 … N) is just the product of the corresponding individual likelihoods: (28) Eq 28 can now be applied within the context of Bayesian data analysis [35, 36, 63] to estimate parameters of Gaussian stochastic models defined by their first moment and covariance in the light of the available experimental data. According to the Bayesian theorem, the posterior distribution of parameters p(θ|d) is given as (29) represents the likelihood as defined in Eq 28. The so-called prior distribution p(θ) contains all information about the model parameters without any relation to the current experimental data set. Z is the so-called evidence: It normalizes the posterior distribution and gives information about the model probability. We apply the nested sampling algorithm [64, 65] as a numerical technique to obtain the marginalized posterior distributions of parameters and model probabilities. Mean parameter values and their uncertainties are contained in the posterior of parameters p(θ|d). They are obtained formally from the calculation of moment integrals over the posterior with respect to the parameters θ1θM (for M parameters): (30) The mean 〈θk〉 of the parameter θk is obtained for n = 1, its standard deviation with n = 2 by (31) These values are given as mean values and error bars for all parameters in Figs 57. Thereby, we use weak constant priors for H ∈ [0.01, 0.99], vd ∈ [−25, 25] μm/min, τ* ∈ [0, 30] min, μ ∈ [0, 4], and an exponential prior ∼ exp(−DH/Dm) for DH with Dm = 50 μm2/min2H [35, 36]. In this work we apply Eqs 28 and 29 to our experimental paths of neutrophils using the stochastic Gaussian models defined by Eqs 810.

The nested sampling algorithm, which is used to calculate the integrals in Eq 30, delivers numerical sampling curves of parameters and evidences until convergence [64, 65]. Each sampling point i (i = 1 … NS with NS being the number of sampling points) receives a normalized weight wi (). These weights are used to calculate mean values and uncertainties of a parameter θk (32) (33) where Eq 31 applies analogously to calculate δθk. In addition, this procedure is also applied to calculate the mean and standard deviations of functions f(θ) depending on the parameters θ: (34) This delivers the real mean value and standard deviation of the function f, which may be different from just calculating the function f for the mean parameters θ. This calculation is applied to obtain the mean logarithmic derivation β(t) of the msd in Figs 5 and 7.

Model probabilities are obtained from the corresponding evidences Z of all considered models A–D. E.g. the probability P(A) of model A is given by (35) Thus, individual model probabilities are between 0 … 1. The definition shows that the sum of the model probabilities of all considered models sums up to one. The resulting model probabilities are shown in Figs 4 and 7 as percentages between 0 … 100%.

Simulations based on the covariance function

Gaussian stochastic processes are completely specified by their mean and covariance. Thus, mean and covariance can be used for simulating trajectories of the corresponding stochastic process. Covariance matrices as given in Eqs 810 are real, symmetric and have positive eigenvalues. Thus, the covariance can be factorized by a Cholesky decomposition [66] into a product of a lower triangular matrix with its transpose . With uncorrelated Gaussian random variables (i = 1…N), N positions of the corresponding trajectory are given as a matrix product with x = {x1, x2, …xN}. This is obvious [67, 68] from the calculation of expectation values as (36) where with the identity matrix has been used. Thus, after performing the Cholesky decomposition of the covariance , sample trajectories of the stochastic process can be generated according to . This method can quite simply be applied when covariances are known [69]. As the numerical efficiency of this procedure is low, specialized algorithms have been developed for specific stochastic processes as fractional Brownian motion (see e.g. [67, 69]). However, for purposes of this work the covariance algorithm was sufficient to perform simulations of the suggested anomalous processes.

The python programs for Bayesian analysis including covariance can be found on the zenodo archive (DOI:10.5281/zenodo.6325568).

Supporting information

S3 Data. Excel file with data shown in Fig 3 and additionally including the velocity cross correlations 〈vcross,xy(t)〉.

https://doi.org/10.1371/journal.pcbi.1010089.s003

(XLSX)

S4 Data. Excel file with data shown in Fig 4 and additionally including the model probabilities of the Ornstein-Uhlenbeck process (model D).

https://doi.org/10.1371/journal.pcbi.1010089.s004

(XLSX)

Acknowledgments

We thank Andreas Deussen and Melanie Martin for discussions and critically reading of the manuscript.

References

  1. 1. Nourshargh S, Hordijk PL, Sixt M. Breaching multiple barriers: leukocyte motility through venular walls and the interstitium. Nat Rev Mol Cell Biol. 2010;11(5):366–378. pmid:20414258
  2. 2. Liew PX, Kubes P. The Neutrophil’s Role During Health and Disease. Physiological Reviews. 2019;99(2):1223–1248. pmid:30758246
  3. 3. Mayadas TN, Cullere X, Lowell CA. The multifaceted functions of neutrophils. Annu Rev Pathol. 2014;9:181–218. pmid:24050624
  4. 4. Cai H, Devreotes PN. Moving in the right direction: how eukaryotic cells migrate along chemical gradients. Semin Cell Dev Biol. 2011;22(8):834–841. pmid:21821139
  5. 5. Steinckwich N, Myers P, Janardhan KS, Flagler ND, King D, Petranka JG, et al. Role of the store-operated calcium entry protein, STIM1, in neutrophil chemotaxis and infiltration into a murine model of psoriasis-inflamed skin. FASEB J. 2015;29(7):3003–3013. pmid:25837581
  6. 6. Beerman RW, Matty MA, Au GG, Looger LL, Choudhury KR, Keller PJ, et al. Direct In Vivo Manipulation and Imaging of Calcium Transients in Neutrophils Identify a Critical Role for Leading-Edge Calcium Flux. Cell Rep. 2015;13(10):2107–2117. pmid:26673320
  7. 7. Dixit N, Simon SI. Chemokines, selectins and intracellular calcium flux: temporal and spatial cues for leukocyte arrest. Front Immunol. 2012;3:188. pmid:22787461
  8. 8. Schaff UY, Yamayoshi I, Tse T, Griffin D, Kibathi L, Simon SI. Calcium flux in neutrophils synchronizes β2 integrin adhesive and signaling events that guide inflammatory recruitment. Ann Biomed Eng. 2008;36(4):632–646. pmid:18278555
  9. 9. Pantaler E, Lückhoff A. Inhibitors of TRP channels reveal stimulus-dependent differential activation of Ca2+ influx pathways in human neutrophil granulocytes. Naunyn Schmiedebergs Arch Pharmacol. 2009;380(6):497–507. pmid:19894037
  10. 10. Damann N, Owsianik G, Li S, Poll C, Nilius B. The calcium-conducting ion channel transient receptor potential canonical 6 is involved in macrophage inflammatory protein-2-induced migration of mouse neutrophils. Acta Physiol (Oxf). 2009;195(1):3–11. pmid:18983454
  11. 11. Lindemann O, Umlauf D, Frank S, Schimmelpfennig S, Bertrand J, Pap T, et al. TRPC6 Regulates CXCR2-Mediated Chemotaxis of Murine Neutrophils. J Immunol. 2013;190(11):5496–5505. pmid:23636057
  12. 12. Najder K, Musset B, Lindemann O, Bulk E, Schwab A, Fels B. The function of TRP channels in neutrophil granulocytes. Pflügers Archiv—European Journal of Physiology. 2018;470(7):1017–1033. pmid:29717355
  13. 13. Iglesias PA, Devreotes PN. Biased excitable networks: how cells direct motion in response to gradients. Curr Opin Cell Biol. 2012;24(2):245–253. pmid:22154943
  14. 14. Pal DS, Li X, Banerjee T, Miao Y, Devreotes PN. The excitable signal transduction networks: movers and shapers of eukaryotic cell migration. The International journal of developmental biology. 2019;63(8-9-10):407–416. pmid:31840779
  15. 15. Neilson MP, Veltman DM, van Haastert PJM, Webb SD, Mackenzie JA, Insall RH. Chemotaxis: a feedback-based computational model robustly predicts multiple aspects of real cell behaviour. PLoS Biol. 2011;9(5):e1000618. pmid:21610858
  16. 16. Harris TH, Banigan EJ, Christian DA, Konradt C, Tait Wojno ED, Norose K, et al. Generalized Lévy walks and the role of chemokines in migration of effector CD8+ T cells. Nature. 2012;486(7404):545–548. pmid:22722867
  17. 17. Amselem G, Theves M, Bae A, Bodenschatz E, Beta C. A stochastic description of Dictyostelium chemotaxis. PLoS One. 2012;7(5):e37213. pmid:22662138
  18. 18. Stokes CL, Lauffenburger DA, Williams SK. Migration of Individual Microvessel Endothelial Cells: Stochastic Model and Parameter Measurement. Journal of Cell Science. 1991;99(2):419–430. pmid:1885678
  19. 19. Codling EA, Plank MJ, Benhamou S. Random Walk Models in Biology. Journal of The Royal Society Interface. 2008;5(25):813–834. pmid:18426776
  20. 20. Codling EA, Bearon RN, Thorn GJ. Diffusion about the Mean Drift Location in a Biased Random Walk. Ecology. 2010;91(10):3106–3113. pmid:21058570
  21. 21. Scott M, Żychaluk K, Bearon R. A Mathematical Framework for Modelling 3D Cell Motility: Applications to Glioblastoma Cell Migration. Mathematical Medicine and Biology: A Journal of the IMA. 2021;38(3):333–354. pmid:34189581
  22. 22. Wu PH, Giri A, Sun SX, Wirtz D. Three-Dimensional Cell Migration Does Not Follow a Random Walk. Proceedings of the National Academy of Sciences. 2014;111(11):3949–3954.
  23. 23. Dieterich P, Klages R, Preuss R, Schwab A. Anomalous dynamics of cell migration. Proc Natl Acad Sci USA. 2008;105(2):459–463. pmid:18182493
  24. 24. Denker SP, Barber DL. Cell migration requires both ion translocation and cytoskeletal anchoring by the Na-H exchanger NHE1. J Cell Biol. 2002;159(6):1087–1096. pmid:12486114
  25. 25. Schwab A, Fabian A, Hanley PJ, Stock C. Role of ion channels and transporters in cell migration. Physiol Rev. 2012;92(4):1865–1913. pmid:23073633
  26. 26. Shi C, Huang CH, Devreotes PN, Iglesias PA. Interaction of motility, directional sensing, and polarity modules recreates the behaviors of chemotaxing cells. PLoS Comput Biol. 2013;9(7):e1003122. pmid:23861660
  27. 27. Tauzin S, Starnes TW, Becker FB, Lam Py, Huttenlocher A. Redox and Src family kinase signaling control leukocyte wound attraction and neutrophil reverse migration. J Cell Biol. 2014;207(5):589–598. pmid:25488917
  28. 28. Andrew N, Insall RH. Chemotaxis in shallow gradients is mediated independently of PtdIns 3-kinase by biased choices between random protrusions. Nat Cell Biol. 2007;9(2):193–200. pmid:17220879
  29. 29. Feller W. Introduction to the Theory of Probability and its Applications, Vol. II. (2 Ed) New York: Wiley. 1971.
  30. 30. Metzler R, Klafter J. The random walk’s guide to anomalous diffusion: a fractional dynamics approach. Physics Reports. 2000;339(1):1–77.
  31. 31. Mandelbrot BB, Ness JWV. Fractional Brownian Motions, Fractional Noises and Applications. SIAM Rev. 1968;10(4):422–437.
  32. 32. Erdélyi A, Magnus W, Oberhettinger F, Tricomi FG. Higher Transcendental Functions (Volume I of the Bateman project). vol. I. New York: McGraw-Hill Book Company; 1953.
  33. 33. Gel’fand I, Shilov G, Saletan E. Generalized functions. Volume 1, Volume 1,. New York: Academic Press Inc.; 1964.
  34. 34. Molina-Garcia D, Sandev T, Safdari H, Pagnini G, Chechkin A, Metzler R. Crossover from anomalous to normal diffusion: truncated power-law noise correlations and applications to dynamics in lipid bilayers. New Journal of Physics. 2018;20(10):103027.
  35. 35. von der Linden W, Dose V, Von Toussaint U. Bayesian probability theory: applications in the physical sciences. Cambridge University Press; 2014.
  36. 36. Gregory PC. Bayesian Logical Data Analysis for the Physical Sciences: A Comparative Approach with Mathematica Support. Cambridge University Press; 2005.
  37. 37. Andreasen E, Mathew L, Löhr C, Hasson R, Tanguay R. Aryl Hydrocarbon Receptor Activation Impairs Extracellular Matrix Remodeling during Zebra Fish Fin Regeneration. Toxicological Sciences. 2007;95(1):215–226. pmid:17003102
  38. 38. Senk A, Djonov V. Collagen Fibers Provide Guidance Cues for Capillary Regrowth during Regenerative Angiogenesis in Zebrafish. Scientific Reports. 2021;11(1):19520. pmid:34593884
  39. 39. Reversat A, Gaertner F, Merrin J, Stopp J, Tasciyan S, Aguilera J, et al. Cellular Locomotion Using Environmental Topography. Nature. 2020;582(7813):582–585. pmid:32581372
  40. 40. Redner S. A Guide to First-Passage Processes. Cambridge University Press; 2001.
  41. 41. Lämmermann T. In the eye of the neutrophil swarm-navigation signals that bring neutrophils together in inflamed and infected tissues. J Leukoc Biol. 2016;100:55–63. pmid:26416718
  42. 42. Lämmermann T, Afonso PV, Angermann BR, Wang JM, Kastenmüller W, Parent CA, et al. Neutrophil swarms require LTB4 and integrins at sites of cell death in vivo. Nature. 2013;498(7454):371–375. pmid:23708969
  43. 43. Schaefer B, Bartosova M, Macher-Goeppinger S, Ujszaszi A, Wallwiener M, Nyarangi-Dix J, et al. Quantitative Histomorphometry of the Healthy Peritoneum. Scientific Reports. 2016;6(1):21344. pmid:26905058
  44. 44. Nielsen N, Lindemann O, Schwab A. TRP channels and STIM/ORAI proteins: sensors and effectors of cancer and stroma cell migration. Br J Pharmacol. 2014;171(24):5524–5540. pmid:24724725
  45. 45. Wei C, Wang X, Chen M, Ouyang K, Song LS, Cheng H. Calcium flickers steer cell migration. Nature. 2009;457(7231):901–905. pmid:19118385
  46. 46. Tsai FC, Meyer T. Ca2+ pulses control local cycles of lamellipodia retraction and adhesion along the front of migrating cells. Curr Biol. 2012;22(9):837–842. pmid:22521790
  47. 47. Tsai FC, Seki A, Yang HW, Hayer A, Carrasco S, Malmersjö S, et al. A polarized Ca2+, diacylglycerol and STIM1 signalling system regulates directed cell migration. Nat Cell Biol. 2014;16(2):133–144. pmid:24463606
  48. 48. Evans JH, Falke JJ. Ca2+ influx is an essential component of the positive-feedback loop that maintains leading-edge structure and activity in macrophages. Proc Natl Acad Sci USA. 2007;104(41):16176–16181. pmid:17911247
  49. 49. Nosbisch JL, Rahman A, Mohan K, Elston TC, Bear JE, Haugh JM. Mechanistic models of PLC/PKC signaling implicate phosphatidic acid as a key amplifier of chemotactic gradient sensing. PLOS Computational Biology. 2020;16(4):1–21. pmid:32255775
  50. 50. Dietrich A, Gudermann T. In: Nilius B, Flockerzi V, editors. TRPC6: Physiological Function and Pathophysiological Relevance. vol. 222. Berlin, Heidelberg: Springer Berlin Heidelberg; 2014. p. 157–188.
  51. 51. Hofmann T, Obukhov AG, Schaefer M, Harteneck C, Gudermann T, Schultz G. Direct activation of human TRPC6 and TRPC3 channels by diacylglycerol. Nature. 1999;397(6716):259–263. pmid:9930701
  52. 52. Nywening TM, Belt BA, Cullinan DR, Panni RZ, Han BJ, Sanford DE, et al. Targeting both tumour-associated CXCR2+ neutrophils and CCR2+ macrophages disrupts myeloid recruitment and improves chemotherapeutic responses in pancreatic ductal adenocarcinoma. Gut. 2017; p. 2017–313738. pmid:29196437
  53. 53. Sharma B, Nannuru KC, Varney ML, Singh RK. Host CXCR2-dependent regulation of mammary tumor growth and metastasis. Clinical & experimental metastasis. 2015;32(1):65–72. pmid:25511644
  54. 54. Tazzyman S, Barry ST, Ashton S, Wood P, Blakey D, Lewis CE, et al. Inhibition of neutrophil infiltration into A549 lung tumors in vitro and in vivo using a CXCR2-specific antagonist is associated with reduced tumor growth. International Journal of Cancer. 2011;129(4):847–858. pmid:21328342
  55. 55. Krummel MF, Bartumeus F, Gerard A. T cell migration, search strategies and mechanisms. Nat Rev Immunol. 2016;16(3):193–201. pmid:26852928
  56. 56. Guérin T, Levernier N, Bénichou O, Voituriez R. Mean first-passage times of non-Markovian random walkers in confinement. Nature. 2016;534(7607):356–359. pmid:27306185
  57. 57. Palyulin VV, Chechkin AV, Metzler R. Lévy flights do not always optimize random blind search for sparse targets. Proc Natl Acad Sci USA. 2014;111(8):2931–2936.
  58. 58. Hou R, Cherstvy A, Metzler R, Akimoto T. Biased Continuous-Time Random Walks for Ordinary and Equilibrium Cases: Facilitation of Diffusion, Ergodicity Breaking and Ageing. Physical Chemistry Chemical Physics. 2018;20(32):20827–20848. pmid:30066003
  59. 59. Mardoukhi Y, Chechkin A, Metzler R. Spurious Ergodicity Breaking in Normal and Fractional Ornstein–Uhlenbeck Process. New Journal of Physics. 2020;22(7):073012.
  60. 60. Makarava N, Benmehdi S, Holschneider M. Bayesian estimation of self-similarity exponent. Phys Rev E. 2011;84:021109. pmid:21928951
  61. 61. Makarava N, Holschneider M. Estimation of the Hurst exponent from noisy data: a Bayesian approach. The European Physical Journal B. 2012;85(8).
  62. 62. Makarava N, Menz S, Theves M, Huisinga W, Beta C, Holschneider M. Quantifying the degree of persistence in random amoeboid motion based on the Hurst exponent of fractional Brownian motion. Phys Rev E. 2014;90:042703. pmid:25375519
  63. 63. Dieterich P, Preuss R. Bayesian parameter estimation for stochastic models of biological cell migration. In: Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 32nd International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering. American Institute of Physics; 2013. p. 16–22.
  64. 64. Skilling J. Nested Sampling for General Bayesian Computation. Bayesian Analysis—International Society for Bayesian Analysis. 2006;1(4):833–860.
  65. 65. Buchner J. A statistical test for Nested Sampling algorithms. Statistics and Computing. 2016;26(1):383–392.
  66. 66. Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical recipes in C: the art of scientific computing. 2nd ed. ISBN 0-521-43108-5. Cambridge University Press; 1992.
  67. 67. Asmussen S, Glynn PW. Stochastic Simulation: Algorithms and Analysis. No. 57 in Stochastic Modelling and Applied Probability. Springer-Verlag New York; 2007.
  68. 68. Dietrich CR, Newsam GN. Fast and exact simulation of stationary Gaussian processes through circulant embedding of the covariance matrix. SIAM Journal on Scientific Computing. 1997;18(4):1088–1107.
  69. 69. Dieker T. Simulation of fractional Brownian motion. University of Twente, Department of Mathematical Sciences, Enschede, The Netherlands; 2004.