Anomalous diffusion and asymmetric tempering memory in neutrophil chemotaxis

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.


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 Ca 2+ concentration by Ca 2+ release from cytosolic stores and Ca 2+ influx via Ca 2+ permeable ion channels in the plasma membrane [5][6][7][8] such as members of the family of transient receptor potential (TRP) channels ( [9][10][11] and [12] for a review). Impairing Ca 2+ 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 [18][19][20][21] 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 [23][24][25]. 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. 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).

Overview of chemotaxis experiments
The red ellipses represent the mean squared propagation of the cell positions in x-and ydirection 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 16-19 (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 v d given as slope of hy(t)i. 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.

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 hv corr,x|y (t)i 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 powerlaw. The velocity autocorrelation contains the diffusive and drift components of the stochastic process: Starting from hv(t = 0) 2 i at time t = 0 it decays towards the mean squared drift velocity v 2 d for long times. In addition, we analyzed the velocity cross-correlations hv cross,xy (t)i according to Eqs 23 and 24 and did not find couplings of velocities between x-and y-direction (see S3 Data).

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 v d (CXCL1 gradient) and an active stochastic force z x|y (t) generated by the cellular migration apparatus: 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) [18][19][20][21][22]. 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.
with the generalized diffusion coefficient D H 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 τ = |t − t 0 |. 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 1=Gð�Þ � � þ Oð� 2 Þ for � ! 0 [32] and a representation of the delta function dðtÞ ¼ lim Thus, the correlation function reduces to hv 2 (τ)i 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 τ � ! 1.
Model C: Fractional Brownian noise with power-law tempering [34] hv with tempering time τ � and power-law tempering exponent μ � 0. Model C is identical with model A for μ = 0 or τ � ! 1.

Model D:
Ornstein-Uhlenbeck like processes are often used as reference processes to describe and quantify cell migration [18][19][20][21][22]. Thus, we also included this approach, which is characterized by an exponentially decreasing velocity autocorrelation function: Parameter v 2 m 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: The second moments given by the covariances of these processes can be obtained by integrating the velocity autocorrelations of Eqs 2-5 twice: They are different for each model: 1 F 1 and 2 F 1 represent the confluent hypergeometric and the Gaussian hypergeometric function, respectively [32]. The y-direction contains additional terms v 2 d t t 0 due to the drift velocity v d , where model parameters (H, D H , τ � , μ) and ðg; v 2 m Þ can take different values with respect to x-and y-direction. The corresponding second moments follow directly from the covariances of Eqs 8-11 for t = t 0 as Analogously to the covariances, the msd in y-direction msd y (t) A|B|C|D is extended by the drift factor v 2 d t 2 . In addition, the mean squared displacements for models A to C in Eqs 12-14 reduce to normal diffusion 2 D 1/2 t for H ! 1/2 as 1 F 1 (0, 2; z) = 1 and 2 F 1 (μ, 0, 2; z) = 1 [32].
The msd x (t) D of the Ornstein-Uhlenbeck process in Eq 15 displays a ballistic phase � v 2 m t 2 for small times t � 1/γ and normal diffusion � 2 v 2 m =g t 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    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).
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), D H (diffusion constant) and v d (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 v d 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 D H (Fig 6D) is about 39.01 ± 1.89 μm 2 /min for WT neutrophils. WT+SB and TRPC6KO neutrophils have a higher diffusion coefficient.
In summary, it is notable that parameters H and D H 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 v d .

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

PLOS COMPUTATIONAL BIOLOGY
Anomalous diffusion and asymmetric tempering memory in neutrophil chemotaxis 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 D H 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 v d � 2.5 μm/min, which is reduced to v d � 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 8-10 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 socalled first passage time t fpt describing the duration neutrophils need to pass a target for the first time [40]. Fig 8A shows the probability of first passage times P(t fpt ) 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.

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 Ca 2+ signaling domains as described earlier [6,[45][46][47][48]. 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 v d 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 v d . 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 D H 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,[55][56][57].

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 x k (t) of cell k as a function of time t in the interval [0, T]. These are defined as 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: The definitions for the y-direction apply analogously. The second moment is synonymously called mean squared displacement in the present work, i.e. msd x (t) = hx 2 (t)i and msd y (t) = hy 2 (t)i. 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 v x|y (t) was calculated from trajectories as the difference quotient of two cell positions at times t + Δ and t.
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 16-19 velocity autocorrelations are defined as v corr;x;k ðtÞ ¼ 1 These definitions apply analogously to the y-direction. Accordingly, the time-and ensembleaveraged velocity cross-correlation function between velocities in x-and y-direction can be defined as v cross;xy;k ðtÞ ¼ 1

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: The mean drift μ and the covariance C may depend on model parameters θ with (compare with Eqs 6-10). The likelihood of a cell ensemble consisting of N paths with data d k (k = 1 . . . N) is just the product of the corresponding individual likelihoods: 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 LðdjCðθÞ 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): The mean hθ k i of the parameter θ k is obtained for n = 1, its standard deviation with n = 2 by  [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 8-10. 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 . . . N S with N S being the number of sampling points) receives a normalized weight w i ( P N S i¼1 w i ¼ 1). These weights are used to calculate mean values and uncertainties of a parameter θ k 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 θ: 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

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 8-10 are real, symmetric and have positive eigenvalues. Thus, the covariance Cðt; t 0 Þ ¼ hxðtÞxðt 0 Þi can be factorized by a Cholesky decomposition [66] Cðt; t 0 Þ ¼ A A T into a product of a lower triangular matrix A with its transpose A T . With uncorrelated Gaussian random variables Z i � N ð0; 1Þ (i = 1. . .N), N positions of the corresponding trajectory are given as a matrix product x ¼ A η with x = {x 1 , x 2 , . . .x N }. This is obvious [67,68] from the calculation of expectation values as where hη η T i ¼ I with the identity matrix I has been used. Thus, after performing the Cholesky decomposition of the covariance C, sample trajectories of the stochastic process can be generated according to x ¼ A η. 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).