Modelling of Dictyostelium Discoideum Movement in Linear Gradient of Chemoattractant

Chemotaxis is a ubiquitous biological phenomenon in which cells detect a spatial gradient of chemoattractant, and then move towards the source. Here we present a position-dependent advection-diffusion model that quantitatively describes the statistical features of the chemotactic motion of the social amoeba {\it Dictyostelium discoideum} in a linear gradient of cAMP (cyclic adenosine monophosphate). We fit the model to experimental trajectories that are recorded in a microfluidic setup with stationary cAMP gradients and extract the diffusion and drift coefficients in the gradient direction. Our analysis shows that for the majority of gradients, both coefficients decrease in time and become negative as the cells crawl up the gradient. The extracted model parameters also show that besides the expected drift in the direction of chemoattractant gradient, we observe a nonlinear dependency of the corresponding variance in time, which can be explained by the model. Furthermore, the results of the model show that the non-linear term in the mean squared displacement of the cell trajectories can dominate the linear term on large time scales.


I. INTRODUCTION
Dictyostelium discoideum (D.d.) is a well-established model organism for cellular motility. Chemotactic competent D.d. cells are highly motile and exhibit fast amoeboid movements with a velocity of 10 − 20 µm/min on glass substrates [1][2][3]. The chemotactic cell motion is highly organized over a length scale significantly larger than the size of a single cell (∼10 µm). When nutrients are depleted, D.d. cells secret a chemical called cAMP (cyclic adenosine monophosphate) that has an attractive effect on the cells themselves. Cells sense gradients of cAMP and direct their chemotactic movements towards regions of higher concentration of cAMP [4]. When chemotactic attraction prevails over diffusion, the chemotaxis can trigger a self-accelerating process until aggregation takes place. As a result, 10 5 -10 6 cells stream towards the aggregation centers and eventually transform into millimeter long slugs and ultimately form fruiting bodies bearing spores for long-term survival and longrange dispersal [5].
Different mathematical models incorporate chemotaxis in different ways; however, a common mechanism is to assume that chemotaxis biases the otherwise random motion of crawling cells along the concentration gradients of chemoattractants [6]. The random cell movement is commonly described as a diffusion and the directional movement along the chemical gradient is incorporated as a combination of diffusion and advection. In the simplest model, the diffusion coefficient and the drift velocity of the cells are assumed to be constant. However, in general, these coefficients depend on both the absolute concentration and the gradient of the chemical [7][8][9][10]. The advection-diffusion equation has been previously used to describe the aggregation phase of D.d. cells where the chemotactic force pulls the amoebas towards the aggregation centers. For example, a model of slime mold aggregation has been introduced by Patlak [6] and Keller [10] in the form of two coupled differential equations. The first equation is an advection-diffusion equation describing the evolution of the concentration of amobae and the second equation is a diffusion equation with terms of source and degradation describing the evolution of the concentration of the signaling molecule. The original form of the Keller-Segel model, would allow the diffusion coefficient and the drift velocity to depend on the cAMP concentration and on the concentration of the amoeba. The case that these coefficients depend on the chemical cencentration but not on the cell density has been considered by Othmer and Stevens [11]. This leads to ordinary mean field Fokker-Planck equations for cell density with space and time dependent coefficients [12]. On the other hand, if we assume that the diffusion coefficient and the drift velocity of the cells depend on their concentration and on the concentration of the secreted chemical, the original Keller-Segel model takes the form of a generalized mean field Fokker-Planck equation.
The statistical characteristics of trajectories of motile D.d. cells have been the subject of several recent studies. These include experiments to characterize chemotactic cell movement in homogenous and inhomogenous chemical cues [7,8,[13][14][15][16][17][18], and parallel theoretical modeling to reproduce statistical features of the experimental observations [7,[16][17][18][19]. Recently, Li et al. have presented an experimental study of the individual cells in a homogeneous medium [18]. They have proposed a generalized Langevin equation for the velocity of individuals. Their data-driven modeling showed a "programmed" periodic motion around a persistent direction of motion on short time scales and ordinary diffusive behavior on long time scales. Moreover, it is also well known that a cAMP gradient induces a bias of the position where pseudopodia emerge [15]. The measured probabilities of pseudopod directions were used to obtain an analytical model for chemotaxis of cell populations [16]. The prediction of the model are similar to measured chemotactic index of wild-type cells as well as the mutants. Besides, although it is well-known that the directed movement of the D.d. cells in response to the chemoattractant cAMP depends both on the absolute value of the local concentration (chemokinesis) and its gradient (chemotaxis), the exact dependency is not well understood.
In this study, we aim to extract the concentration dependencies of the diffusion and drift coefficients in the Fokker-Planck equation (with respect to cell density), by analyzing the experimental trajectories of motile D.d. cells in Ref. [7,8,13]. We assume that these coefficients depend on both the local cAMP concentration (the so-called midpoint concentration) and its gradient. The experiments are performed in a microfluidic device (see Section II-B) that generate linear stable gradients between the two inlet concentrations C max and C min . As the cells crawl up the gradient, the average background concentration they experience increases. These experiments systematically explore different steepnesses and cover a wide range of gradients, at which chemotactic behavior is observed. In these experiments, an external flow removes the naturally produced cAMP secreted by the cells to avoid cell-cell signaling. This is completely different from aggregation process where the cell density is much higher and the cells signal each other. We start our analysis by assuming linear dependencies for the diffusion coefficient and the drift velocity of the cells along the width of microfluidic setup, where a linear gradient is established. We then use the experimental cell trajectories to deduce the coefficients of these linear dependencies at different cAMP gradients.

A. Cell Culture
All experiments were performed by M. Theves [7,8,13] with Dictyostelium discoideum AX3 wild type cells. Cells were grown in HL5 medium (7 g/L yeast extract, 14 g/L peptone, 0.5 g/L potassium dihydrogen phosphate, 0.5 g/L disodium hydrogen phosphate, 13.5 g/L glucose, ForMedium Ltd., UK). Cells were starved in shaking suspension of phosphate buffer (pH 6.0, 15 mM KH2 PO4 , 1 mM Na2 HPO4 ) at a density of 2 × 10 6 cells/mL for 5:30 hours. After one hour of starvation, the cells were exposed to periodic pulses of cAMP for the remaining time of starvation. The pulses had a concentration of 50 nM and were delivered with a period of 6 minutes.

B. Microfluidics
A microfluidic gradient mixer [17,20] with given dimensions (width=525 µm, height= 50 µm) was used to establish a stable linear gradient over a region of 350 µm × 50 µm × 3000 µm in size (see Fig. 1). The gradients were generated using a pyramidal microfluidic network that provides well-defined concentration profiles with high temporal stability. Throughout the experiment, a constant flow is provided by syringe pumps. The flow provides a constant supply of oxygen and removes all substances released by the cells. This prevents cells from signaling each other, which would perturb the concentration gradient and bias the chemotactic motion. Running at an adjustable average flow velocity ofv = 320 µm/s, the gradient is linear and stable within d = 350 µm in the middle of the channel. Above a lower threshold of ∇C thresh ∼ 10 −3 nM/µm cells started to show a directional response. It is important to note, that all gradients have been established by mixing a phosphate buffer solution at one inlet, C min = 0, together with a solution of cAMP and phosphate buffer C max on the opposing inlet. Therefore the gradient always ranges from zero to this maximum concentration. Due to boundary effects, the profile is distorted near the walls. All cell trajectories within this non-linear area were excluded from statistics. Moreover, given the dimension of the channel and the dynamics viscosity of the flowing phosphate buffer (η = 10 −3 Pa s), one can calculate the shear stress applied on the cells at the imposed mean flow velocity ofv = 320 µm/s to be σ = 0.038 Pa. According to the literature, mechanosensing in D.d. has been observed above a threshold of σ = 0.5 Pa [21].
We are thus approximately one order of magnitude below the regime where flow induced shear stress would bias the motion of chemotactic cells.

C. Cell Tracking
Differential Interference Contrast (DIC) images were recorded for 180 min, with time resolution of 10 sec and spatial resolution of 1024x1024 pixel (1 pixel=0.6409 µm), and processed using Mathworks MATLAB 7.5 with the Image Processing Toolbox [7,8,13]. All the image processing steps are done by M. Theves et al. The images were binarized to distinguish the cells from the background and possible optical artifacts. The cell centroids in each binarized frame were identified. To produce cell trajectories one had to link these locations together in time and space. To achieve this, a customized version of the MATLAB cell tracking algorithm written by Crocker and Grier [22] was used. This tracking pro- cess is consisted of calculating and minimizing the sum over the squared displacements of all possible links between the cell positions in two subsequent frames. Interestingly, an analysis of the broken trajectories have shown that more than 90 percent of the cells were lost due to a sudden jump in the cell location or because two cells ran into each other and their center of mass in the binarized image became indistinguishable. In this case, the tracks will end and new ones will start, once the cells separate again. Tracks may also end when the segmentation program loses a cell due to image quality problems. Once the cell is detected again, a new track will start. These different scenarios result in a distribution of tracks of different length with most of them shorter than the total measurement time. They also result that the number of trajectories (e.g. 582 trajectories in Fig. 2) is much greater than the number of cells (∼ 40 cells at t = 0) in the experiment. Note that during an experi-mental recording, the number of cells is not constant with time as most of the fast cells move out of the region of interest or new cells enter the field of view. Finally, it is important to note that since the cells begin responding to the cAMP at different time points, or as the cells collide and new tracks start, the starting time-points of all trajectories are set to zero.

D. Selection of Trajectories
In our analysis, to have a reliable statistics, we keep the number of trajectories during the averaging process constant. Trajectories are selected based on two criteria: (i) they should persist at least 20 min and (ii) within this time interval, the cells should migrate more than 20 µm in −ŷ direction. The minimum displacement of 20 µm in the direction of gradient, for t = 20 min, gives an average motility ofv y > 1 µm/min. Cells withv y < 1 µm/min are neglected to exclude dead or immobile cells from statistics. As previously mentioned, we lose track of the cells once they collide. Therefore, it is important to note that based on this criteria, if a cell collides with another cell and the time interval between two successive collisions is less than 20 min, this trajectory is excluded from statistics although the cell was crawling withv y > 1 µm/min. Eventually, to improve our statistics, long trajectories, are truncated at 20, 40, 60,... min and, if the conditions above are satisfied, trajectories between 20 to 40 min, 40 to 60 min,etc. are considered as new trajectories and the starting time point of each trajectory is set to zero (see Fig. 2).

III. MODEL
Nonlinear mean field Fokker-Planck equations can find important applications in the context of chemotaxis [23]. Here, we attempt to implement an advection-diffusion approach to describe the chemotactic movement of the D.d. cells experiencing a linear stationary gradient [7,8]. The statistical properties of the system are characterized by the values of the model parameters returned after fitting the model to the experimental trajectories.
To study the chemotactic movement of the D.d. cells, we consider an advection-diffusion model in which the centroid of the cell's perimeter is represented as the position of a particle. We define an orthonormal basis with the unit vectorsx andŷ, wherex is the flow direction and −ŷ is the direction of the spatial gradient of cAMP (see Fig. 2). The position of each cell is given by r = xx + yŷ. The concentration of the D.d cells is low enough, so we can assume that each cell does not sense the presence of the other cells. As stated in Experiment section, microfluidic gradient along y direction is generated by a continuous flow along x. To avoid mixing up the issues of chemotaxis in response to the chemoattractant and mechanotaxis under the influence of the shear stress due to viscous forces, we limit our model to the chemotactic movement of the cells along y. Let us assume that p(x, y, t) denotes the number density of cells at position (x, y) at time t. Then, we have the probability density P (y, t), which is the original density p(x, y, t) integrated over x. The current density along y, J, reads as where v and D are drift velocity and diffusion coefficient, respectively, and ∂ y means differentiation with respect to y. Now, the continuity equation for P and J then reads as where ∂ t denotes differentiation with respect to t. Using Eqs. (2) and (3), one can find the diffusion-advection equation for the problem as The chemotactic motion of the cells depends on both the absolute local concentration (chemo-kinesis) and its gradients (chemotaxis) [7,17]. Based on the experiments of Ref. [7], here we consider a constant spatial gradient of cAMP in the direction of −ŷ. Therefore, one can expect that both the diffusion coefficient, D, and the drift velocity, v, depend on the y component of the position vector. Since there is no direct experimental method to determine this dependency, it is plausible to expand the mentioned coefficients in terms of y as and We keep the terms only up to the first order of y, and treat v 1 , and D 1 as perturbation coefficients. Here, we assume that the current in the y direction does not depend on x. Using equations (4) to (6), one finds The mean value of y-component of the cells' positions is obtained by Differentiating the above expression with respect to time results After substituting ∂ t P (y, t) from Eq. (7) and integrating, one can find By solving this simple ordinary differential equation we find where y 0 ≡ y | t=0 denotes the mean initial y-position of the cells. As it has been mentioned above, v 1 and D 1 are the small parameters and in our model, they have been considered as perturbation parameters. After expanding the exponential factor and keeping the terms up to the first order of perturbation parameters, v 1 and D 1 , one can find It is worth mentioning that since terms like v 1 D 1 are the second order of perturbation parameters, these terms are dropped. The variance of the cells' positions along y is defined as σ 2 y (t) ≡ y(t) 2 − y(t) 2 . Using similar method (see Appendix I for details), one can find σ 2 y (t) as where σ 2 y (0) is the initial variance of the cells' positions along y. We note that in Eq. (13), we have kept the terms up to the first order of perturbation parameters as well.

IV. RESULTS
Now we are in a position to determine the perturbation parameters of our model based on the experimental trajectories. The mean displacement of chemotactic cells and the corresponding variance can be calculated from the experimental trajectories as defined in Appendix I (Eqs. [23][24][25][26]. To characterize the chemotactic behavior of D.d. cells, based on Eqs. (12) and (13), we need to determine the values of v 0 , v 1 , D 0 and D 1 . We treat these factors as tuning parameters and find their values simultaneously by fitting (MATLAB, R2016b) the model to the experimental values of y exp and σ 2 y,exp . Tables I-III include the best fit values of the above mentioned parameters at different cAMP gradients. In Table I Table. II).
• The mean square displacement function σ 2 y (t) shows decreasing behavior at C max = 10, 50, 316 and 10000 nM . This trend is an experimental observation, independent of the introduced model, and has to do with the fact that the cells tend to migrate to areas with high cAMP concentration (top part of the channel). Since in these areas midpoint cAMP concentration is high and most of the cAMP receptors are saturated, therefore the cells slow down and accumulate at the top of the channel. This "accumulation" can give rise to a decreasing σ 2 y (t). In the other word, the possible decrease in the variance is due to a drift towards the top areas of the channel.
• The diffusion coefficient in y direction D 0 + D 1 y is initially positive for all concentrations but as the cells migrate upwards and the value of y decreases, it becomes negative for all concentrations except for C max = 1 nM . We think that this negative diffusion coefficient extracted from the data is an artifact of the perturbative approximation.
To be more exact, the diffusion coefficient depends on the position, and we have Taylor-expanded it and kept only the zeroth and the first order terms (the latter as a perturbative term). While the full position-dependent diffusion coefficient should probably be non-negative, there is no such restriction on its truncated form (which contains only the first two terms): it is just a parameter which is determined through a best fit. Of course the pa-rameters should respect the non-negativity of the variance, and they do, as it is seen that the fitted variance does not become negative.
• We define the mean drift velocity of chemotactic D.d. cells as which shows that the drift velocity depends not only on v 0 and v 1 but also on D 1 . The coefficient v 0 v 1 is a small number for different cAMP gradients (see Table. II), therefore v drift is essentially constant in time. The extracted values of drift velocity at time zero are listed in the fifth column of Table. II. It is interesting that the drift velocity in y direction doesn't depend significantly on the cAMP gradient and fluctuates around 4 µm/min. This is consistent with an independent data analysis performed by M. Theves [13] (see Fig. 3): within a plateau, ranging from 10 −2 nM/µm ≤ ∇C ≤ 1 nM/µm over two orders of magnitude, the chemotactic velocity is almost constant. For gradients above ∇C = 1 nM/µm the directionality of movement is decreased. Exceeding an upper threshold of ∇C up ∼ 10 2 nM/µm, the cell motion becomes isotropic again.
• For all gradients, while the cells crawl up the gradient, the magnitude of v y = v 0 + v 1 y decreases as the midpoint concentration increases. However, the independent data analysis of M. Theves [13] shows a transition: for shallow gradients, right after the onset of chemotaxis ∇C = 0.003 nM/µm, v y increases as the background concentration rises. This effect reverses for steep gradients ∇C > 0.3 nM/µm (see right panel of Fig. 3).

V. DISCUSSION
We have analyzed large data sets of D.d. chemotaxis in linear gradients of cAMP recorded by Theves et al. in a microfluidic setup [7,8,13]. Data sets with different steepnesses of cAMP gradient were included in our analysis, covering a large range of gradients, in which chemotactic behavior was observed. Inspired by the experimental conditions of Ref. [7], we introduced a minimal phenomenological model that explicitly incorporated the dependency of diffusion matrix and velocity of the cells on their positions which corresponds to the position dependence of the local concentration of chemotactic cues. Based on this model, we extracted the physical properties of the chemotactic D.d. cells using the mean and variance of the experimental cell tracks. What is the benefit of this phenomenological model? As highlighted  [13] showing (left) chemotaxis as a function of gradient steepness ∇C: above a threshold at ∇C 10 −3 nM/µm cells show a positive (in our coordinate system negative) average velocity in gradient direction (vy) as well as an increased total motility v, while the perpendicular velocity component in flow direction (vx) remains random and averages to zero within error bars. For gradients ranging over two orders of magnitude, 10 −2 nM/µm ≤ ∇C ≤ 1 nM/µm, the chemotactic speed is constant. (right) Average chemotactic velocity vy as a function of gradient steepness evaluated separately for three different areas, subdividing the region of interest (see Fig. 1d). The midpoint concentration decreases from bin 1 to bin 3. For shallow gradients, the chemotactic velocity increases with a raise in the average midpoint concentration. This effect seems to reverse for steep gradients above 1 nM/µm, where vy decreases slightly in higher concentration backgrounds. The figures are used by the courtesy of M. Theves from his master thesis [13].  previously, chemotactic movement of the cells depend on both the chemoattractant gradient and the average ambient chemoattractant concentration (midpoint concentration). In the microfluidic setup of Theves et al. [7], the cells are exposed to a constant gradient, while the midpoint concentration increases when the cells are moving up the gradient. Traditionally, chemotactic cell motion is described by Langevin-type equation where for each cell track, the velocity and acceleration of the cells are calculated at each point by finite differences from the cell positions [7,8,24]. Therefore in these types of analysis midpoint concentration is globally averaged out. Other quantities such as chemotactic index, defined as the distance moved in gradient direction divided by the total distance, are also averaged quantities where information on mid-point concentration is lost. However in our analysis, instead of velocity and acceleration, we work directly with spatial position of the cells and explicitly include the dependence of the diffusion coefficient and the drift velocity on the midpoint concentration. Taylor expansion of these coefficients up to the first order in y leads to a closed set of equations that can be solved to obtain the fitting parameters. It is worth to check the effect of the dependency of the diffusion and velocity on the local concentration. To do this, let us assume that the drift velocity and the diffusion coefficient were constant. We denote the constant drift velocity and diffusion constant byṽ 0 andD 0 , respectively, to obtain These equations predict a linear dependency on t, in both the mean and the variance of the position, which is not   consistent with the experimental data especially in the variance of σ 2 y (see Fig. 2 and Figs. 6-11). In fact it has been shown in Ref. [25], that any linear diffusion model (even anomalous), which enjoys both time translational invariance and space translational invariance leads to means and variances which are at most linear in time. This is a motivation to use drift and diffusion parameters which do depend on the position. There is an obvious position-dependence in our system: C (the concentration of cAMP) does depend on the coordinate y. Assuming that the drift and diffusion parameters do depend on C, one is left with a y-dependence in the drift and diffusion parameters. A simple manageable model is to Taylorexpand this y-dependence and keep only terms which are up to first order in y. The result is a first order perturbation model, which has been studied here. We emphasize that the experimental conditions of Ref. [7] fulfill the necessary conditions of the mentioned study.
In previous studies, wild-type and mutated epithelial canine kidney cells, it has been shown that the cell dynamics can be characterized by an anomalous diffusion [26]. In particular, mean squared displacement   x as a function of t, before and after shifting all the tracks to x = 0, respectively. Red lines correspond to the experimental data and blue lines correspond to fitted quadratic polynomials, respectively.
shows a super-diffusive behavior. This super-diffusive behaviour was also observed in the mean square displacement of Hydra cells [27]. However, experimental trajectories of chemotactic D.d. cells in Ref. [18], were interpreted by a data-driven model with purely diffusive behavior. As we mentioned above, a pure diffusive model can not explain the non-linear behavior observed in the experimental data of Ref. [7].
In our analysis, we observed that at all concentrations D 0 +D 1 y decreases with time and becomes negative for concentrations of 10, 31.6, 316, and 10000 nM . In order to make sure that negative diffusion coefficients are not due to our low statistics after the selection procedure, we combined trajectories of three different experiments performed at the same cAMP gradient, namely ∇C = 0.14 nM/µm (C min = 0, C max = 50 nM ). Comparison between Fig. 2 and Fig. 4 shows a similar decreasing behavior in σ 2 y (t). Indeed, an absorbing point on top of the channel can produce a decreasing variance, not through the diffusion but through the upstream drift. Let us suppose that D 0 = D 1 = 0. Then, according to Eq. 21, σ 2 y (t) = σ 2 y (0) exp(2v 1 t). If everything is expanded up to first order in v 1 , then the result is σ 2 y (t) = σ 2 y (0)(1 + 2v 1 t).
As v 1 ≤ 0, it seems that σ 2 y (0) could become negative after a while. But that is an artifact of the approximation. We intended to find positiondependence of diffusion and drift coefficients. Since the exact position-dependence is not known, even if the inhomogeneity of the surface is known, we expanded the diffusion and drift coefficients in power of the position y. That the time dependence of the variance does match the experiments, means that the method works. But the perturbative parameters should not be misleading.
Furthermore, with our model, we can test directly the space-time symmetries of the cell movement. Based on the reports of the experiments the gradient in the y direction is homogeneous in x. However, the cell tracks shown in the panel (a) of Figs. 2,6,7,9,10, and 11 seem to show a drift in positive x direction (in addition to the chemotactic drift in -y direction). To check the spatial homogeneity in the x direction, we shifted all the tracks of Fig. 2 to x = 0 (see Fig. 5). It is interesting that in this case σ 2 x (t) is not a pure translation of the same function for unshifted trajectories (see Fig. 5). It seems that the behavior of the function depends on the initial condition. This non-pure shift in σ 2 x (t) could be a hallmark of correlation between the displacement of individual cells along x direction and their initial xpositions (see Appendix II). Apparently, the system does not have the translational symmetry along the x direction. This is surprising, since analysis in Section II.B. shows that with flow speed of 320 µm/s we are far below the regime, where mechanotactic effects have been observed in D.d. cells. However, the authors of Ref. [21] conducted their experiments with vegetative cells. This suggests that starvation may increase the mechanosensitivity of D.d. cells. We emphasize that with this correction all of our analysis in the y direction is still valid, given that the current in the y direction does not depend on x. This assumption is nothing but a mean-field approximation.
To improve our statistics, we have divided long mother trajectories to shorter ones and if the criteria in Section. II.D are satisfied, we have included daughter trajectories as completely independent tracks in our analysis. The main difference between these new daughter trajectories is the average midpoint concentration that the cells experience as they crawl up the gradient. This corresponds to moving up from Bin 3 to Bin 1 in Fig. 1d, where in each Bin cells are exposed to a different average midpoint concentration. Detailed analysis by Theves et al. have shown that (see the right panel of Fig. 3) with a raise in the average midpoint concentration the average chemotactic velocity v y doesn't show any clear trend for intermediate gradients, 10 −2 nM/µm ≤ ∇C ≤ 0.3 nM/µm. However, for shallow gradients v y increases with midpoint concentration and for steep gradients it decreases. Most of our analysis are done at intermediate gradients where the variation in chemotactic velocity v y between three different Bins is less than 25 percent. Exemplary, at ∇C = 0.3 nM/µm, average chemotactic velocity changes from ∼ 3 µm/min to ∼ 4 µm/min for three bins with different midpoint concentrations of 17 nM, 50 nM and 83 nM . At steep gradients larger than 0.3 nM/µm, the variations in v y is even less than 10 percent. Thereby we believe that shorter daughter trajectories which be-long to one mother long trajectory, do not significantly differ in their chemotactic properties. In other words, by dividing long mother trajectories to shorter daughter tracks, we don't introduce new types of trajectories with completely different statistical properties.
In the present work, even though we have analyzed a substantial amount of data, much larger data sets with longer trajectories would be required in order to improve our statistics. Possible future experiments with wider microfluidic channels can be helpful to obtain longer trajectories. Experiments with lower cell density (to avoid cell-cell collision) can also help us to obtain longer trajectories, as the cells after collision are indistinguishable from each other and two new trajectories are detected by the cell tracking algorithms.
In summary, we have analyzed the experimental data of chemotactic D.d. cells in the linear gradient of cAMP. In order to have a reliable statistics, we kept the number of trajectories during our analysis constant. Trajectories were selected based on two criteria: (i) they should persist at least 20 min, and (ii) within this time interval, the cells should have migrated more than 20 µm in the direction of the the gradient of cAMP. We have shown that by introducing an advection-diffusion model that includes the position dependence on the cAMP concentration, a quantitative description of experimental cell tracks of the amoeba D.d. is achieved. Our analysis goes beyond a pure diffusive model and shows that the super-diffusive behavior can dominate at larger time scales. Specifically, while in a conventional advection-diffusion model both the mean and the variance are linear in time, here in both cases terms arise which are quadratic in time.
In future study, we aim to apply our analysis to the trajectories of cells migrating on surfaces of differing composition. In a recent study, it has been shown that D.d. cells migrate similarly on surfaces with various chemical composition [28]. As the substrate composition changes, the cells regulate forces generated by actomyosin network to maintain optimal cell-surface contact area and adhesion. We will assess migration trajectories of the cells on different surfaces and investigate the variations in the fitting parameters of our model. Furthermore, we aim to extend our analysis to the trajectories of mutant cell lines that single or multiple components of the chemotactic signaling pathway are deficient and consequently the character of the cell trajectories may change considerably. Structural differences between the trajectories of wild-type and mutant cells may reflect important information about the role of the various proteins in the signaling pathway of D.d. cells, which possibly can not be resolved in the models that mid-point concentration informations are averaged out. The objective is to correlate various parameters of our model to the key molecular players involved in chemotaxis. This can provide a link between the observed macroscopic dynamics and the underlying microscopic mechanism which is an important goal in the field of eukaryotic chemotaxis.

Appendix I
Here we show how to derive Eq. (13). Indeed, for what follows we do not need the exact form of P (y, t) itself, but just the time dependence of its moments. y(t) 2 is defined as y(t) 2 = y 2 P (y, t)dy.
We can directly obtain an equation for the time evolution of y(t) 2 by multiplying the master equation, Eq. (7), by y 2 and integrate over y. This results in ∂ ∂t dy y 2 P (y, t) = dy y 2 (D 0 + D 1 y) ∂ 2 y P The left-hand side of this equation is simply equal to d y(t) 2 dt . We apply partial integration to the right-hand side and obtain d dt y(t) 2 = 2D 0 + 2(v 0 + 2D 1 ) y(t) + 2v 1 y(t) 2 .
Since σ 2 y (t) ≡ y(t) 2 − y(t) 2 , and d dt σ 2 y (t) = d dt y(t) 2 − 2 y(t) d dt y(t) , one finds d dt σ 2 y (t) = 2D 0 + 2D 1 y(t) + 2v 1 σ 2 y (t), where according to Eq. (10), d y /dt has been replaced by v 0 + D 1 + v 1 y(t) . The solution of Eq. (20) is found as σ 2 y (t) = σ 2 y (0) + where σ 2 y (0) denotes the initial variance of the cells' positions along y. After expanding the above equation and keeping the terms up to the first order of v 1 and D 1 , one has σ 2 y (t) = σ 2 y (0) + 2 σ 2 The mean displacement of chemotactic cells, in both x and y directions, and the corresponding variances can be calculated from the experimental data as follows where N denotes the number of cells.