DOI:
10.1039/C7SM01568B
(Paper)
Soft Matter, 2017,
13, 82098222
Modelling of Dictyostelium discoideum movement in a linear gradient of chemoattractant
Received
4th August 2017
, Accepted 27th September 2017
First published on 2nd October 2017
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 positiondependent advection–diffusion model that quantitatively describes the statistical features of the chemotactic motion of the social amoeba 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 over 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 the chemoattractant gradient, we observe a nonlinear dependency of the corresponding variance on time, which can be explained by the model. Furthermore, the results of the model show that the nonlinear 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 wellestablished 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^{−1} on glass substrates.^{1–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 secrete 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 selfaccelerating 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 longterm 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–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 homogeneous and inhomogeneous chemical cues,^{7,8,13–18} and parallel theoretical modeling to reproduce statistical features of the experimental observations.^{7,16–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 datadriven 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 is similar to the measured chemotactic index of wildtype cells as well as the mutants. Besides, although it is wellknown 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 and 13. We assume that these coefficients depend on both the local cAMP concentration (the socalled midpoint concentration) and its gradient. The experiments are performed in a microfluidic device (see Section II.B) that generates 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 an 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 the 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.
II. Experiments
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^{−1} yeast extract, 14 g L^{−1} peptone, 0.5 g L^{−1} potassium dihydrogen phosphate, 0.5 g L^{−1} disodium hydrogen phosphate, 13.5 g L^{−1} glucose, ForMedium Ltd, UK). Cells were starved in shaking suspension of phosphate buffer (pH 6.0, 15 mM KH_{2}PO_{4}, 1 mM Na_{2}HPO_{4}) at a density of 2 × 10^{6} cells per 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 welldefined 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 of = 320 μm s^{−1}, 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^{−1} 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 
∇C = C_{max} − C_{min} = ΔC/d  (1) 
always ranges from zero to this maximum concentration. Due to boundary effects, the profile is distorted near the walls. All cell trajectories within this nonlinear area were excluded from statistics. Moreover, given the dimension of the channel and the dynamic 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 of = 320 μm s^{−1} to be σ = 0.038 Pa. According to the literature, mechanosensing in D.d. cells 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.

 Fig. 1 (a) Microfludic gradient mixer for chemotaxis experiments: two different concentrations flowed into the channel inlets undergo steps of diffusive mixing at each branch to form a linear stable gradient in the area of observation. (b and c) Line profile of fluorescein intensity inside the gradient mixer perpendicular to the flow direction, (d) differential interference contrast (DIC) image, showing the cell population being exposed to the gradient. Only cell trajectories within the region of interest (blue box) are considered for statistics. Moreover, bin 1 corresponds to the area with the highest, and bin 3 to the one with the lowest average midpoint concentration experienced by the cells. This figure is used by the permission of M. Theves from his master's thesis.^{13}  
C. Cell tracking
Differential interference contrast (DIC) images were recorded for 180 min, with time resolution of 10 s and spatial resolution of 1024 × 1024 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 were carried out 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 process 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 has 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 in the number of trajectories (e.g. 582 trajectories in Fig. 2) being much greater than the number of cells (∼40 cells at t = 0) in the experiment. Note that during an experimental 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 timepoints of all trajectories are set to zero.

 Fig. 2 (a) 582 trajectories tracked in a microfluidic channel with C_{max} = 50 nM (∇C = 0.14 nM μm^{−1}). Cells migrate on average upwards from the bottom of the channel to the top areas with higher cAMP concentration. (b) 88 trajectories selected (out of 582) based on the two conditions explained in Section II.D. The stars mark the cell positions exactly at 20 min and (if the trajectory is long enough) at 40 min, 60 min etc. (c) The same trajectories in (b) truncated at 20 min to keep the number of cells during the averaging process constant. For long trajectories, if the two conditions are satisfied, the tracks between 20 to 40 min or from 40 to 60 min, etc., are considered as new independent trajectories to improve the statistics. (d and e) The comparison between experimental data (red lines) and the fitted model (blue line) for 〈y〉 and σ_{y}^{2}.  
D. Selection of trajectories
In our analysis, to have reliable statistics, we keep the number of trajectories during the averaging process constant. Trajectories are selected based on two criteria: (i) they should persist for at least 20 min and (ii) within this time interval, the cells should migrate more than 20 μm in the −ŷ direction. The minimum displacement of 20 μm in the direction of the gradient, for t = 20 min, gives an average motility of _{y} > 1 μm min^{−1}. Cells with _{y} < 1 μm min^{−1} are neglected to exclude dead or immobile cells from the 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 with _{y} > 1 μm min^{−1}. Eventually, to improve our statistics, long trajectories, are truncated at 20, 40, 60,… min and, if the conditions above are satisfied, trajectories between 20 and 40 min, 40 and 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 vectors and ŷ, where 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 = x + 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 the Experimental section, a microfluidic gradient along the 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 the 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
eqn (2) and (3), one can find the diffusion–advection equation for the problem as

∂_{t}P = ∂_{y}(D∂_{y}P) − ∂_{y}(vP).  (4) 
The chemotactic motion of the cells depends on both the absolute local concentration (chemokinesis) 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

v = v_{0} + v_{1}y + …  (5) 
and

D = D_{0} + D_{1}y + …  (6) 
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
eqn (4)–(6), one finds

∂_{t}P = (D_{0} + D_{1}y)∂_{y}^{2}P + (D_{1} − v_{0} − v_{1}y)∂_{y}P − v_{1}P.  (7) 
The mean value of the ycomponent of the cells' positions is obtained by

 (8) 
Differentiating the above expression with respect to time results

 (9) 
After substituting ∂_{t}P(y,t) from eqn (7) and integrating, one can find

 (10) 
By solving this simple ordinary differential equation we find

 (11) 
where 〈
y〉
_{0} ≡ 〈
y〉
_{t=0} denotes the mean initial
yposition of the cells. As it has been mentioned above,
v_{1} and
D_{1} are 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

 (12) 
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 σ_{y}^{2}(t) ≡ 〈y(t)^{2}〉 − 〈y(t)〉^{2}. Using a similar method (see Appendix I for details), one can find σ_{y}^{2}(t) as

σ_{y}^{2}(t) = σ_{y}^{2}(0) + 2[σ_{y}^{2}(0)v_{1} + D_{0} + D_{1}〈y〉_{0}]t + (2D_{0}v_{1} + D_{1}v_{0})t^{2},  (13) 
where
σ_{y}^{2}(0) is the initial variance of the cells' positions along
y. Please note that in
eqn (13), we have kept terms up to the first order of perturbation parameters.
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 (eqn (23)–(26)). To characterize the chemotactic behavior of D.d. cells, based on eqn (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 1–3 include the best fit values of the above mentioned parameters at different cAMP gradients. In Table 1, 〈y〉_{0} and σ_{y}^{2}(0) denote the fitted mean and standard deviations at time zero.
Table 1 The mean initial positions of the cells and the corresponding variances at time zero for different cAMP concentrations
C
_{max} (nM) 
∇C (nM μm^{−1}) 
〈y〉_{0} (μm) 
σ
_{
x
}
^{2}(0) (μm^{2}) 
σ
_{
y
}
^{2}(0) (μm^{2}) 
1 
0.003 
373.39 
10880 
5454.2 
10 
0.03 
347.22 
25893 
6857.1 
31.6 
0.09 
417.88 
24136 
7608.9 
50 
0.14 
365.99 
27253 
6707.3 
100 
0.29 
410.24 
32453 
4719.2 
316 
0.9 
371.13 
30854 
7033.1 
10000 
28.6 
382.33 
27665 
6522.6 
Table 2 Drift coefficients for different cAMP gradients. The mean drift velocity of the cells at time zero is presented in the fifth column
C
_{max} (nM) 
v
_{0} (μm min^{−1}) 
v
_{1} (min^{−1}) 
v
_{0} + v_{1}〈y〉_{0} 
v
_{0} + v_{1}〈y〉_{0} + D_{1} 
v
_{0}
v
_{1}/2 
1 
2.87 
−0.017 
−3.32 
−2.94 
−0.024 
10 
−7.62 
−0.012 
−11.69 
−4.06 
0.045 
31.6 
−7.27 
−0.010 
−11.04 
−4.66 
0.033 
50 
2.04 
−0.016 
−3.95 
−3.45 
−0.017 
100 
3.36 
−0.016 
−3.16 
−3.08 
−0.027 
316 
−8.98 
−0.013 
−13.77 
−4.83 
0.058 
10000 
−10.71 
−0.015 
−16.36 
−3.97 
0.079 
Table 3 Diffusion coefficients in the y direction for different cAMP gradients
C
_{max} (nM) 
D
_{0}

D
_{1}

D
_{0} + D_{1}〈y〉_{0} 
2D_{0}v_{1} + D_{1}v_{0} 
2(σ_{y}(0)^{2}v_{1} + D_{0} + D_{1}〈y〉_{0}) 
1 
−49.10 
0.37 
89.50 
2.69 
−1.81 
10 
−2612.60 
7.63 
38.23 
3.09 
−84.38 
31.6 
−2581.40 
6.38 
86.67 
0.27 
35.76 
50 
−148.60 
0.50 
32.52 
5.87 
−154.38 
100 
129.80 
0.09 
164.54 
−3.84 
179.00 
316 
−3270.57 
8.90 
45.22 
4.15 
−91.00 
10000 
−4690.80 
12.39 
46.59 
6.00 
−99.65 
Fig. 2 and 6–11 (see Appendix III) show the comparison between the model and the experiments at different cAMP concentrations. The initial number of trajectories before the selection procedure is presented in part (a) of each figure. Trajectories for our analysis are then selected and truncated based on the criteria explained in Section II.D. Selected and truncated trajectories are shown in parts (b and c) of each figure, respectively. The initial number of trajectories as well as the number of selected trajectories are different for different cAMP gradients. In parts (d and e) of each figure, the red lines correspond to the experimental data and the blue lines correspond to the results of our model using the fitting parameters of Tables 1–3. The important features of the figures and the tables are summarized below:
• The mean position of the chemotactic D.d. cells, 〈y〉, decreases almost linearly over time, which shows that the chemotactic cells migrate towards areas with higher cAMP concentration (top areas of the channel). The nonlinear term v_{0}v_{1}/2 in eqn (12) is small for all concentrations (see the last column of Table 2).
• The mean square displacement function σ_{y}^{2}(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 the midpoint cAMP concentration is high and most of the cAMP receptors are saturated, the cells slow down and accumulate at the top of the channel. This “accumulation” can result in a decreasing σ_{y}^{2}(t). In other words, the possible decrease in the variance is due to a drift towards the top areas of the channel.
• The diffusion coefficient in the 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 Taylorexpanded it and kept only the zeroth and the first order terms (the latter as a perturbative term). While the full positiondependent diffusion coefficient should probably be nonnegative, 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 parameters should respect the nonnegativity 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

 (14) 
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 2), therefore
v_{drift} is essentially constant over time. The extracted values of drift velocity at time zero are listed in the fifth column of
Table 2. It is interesting that the drift velocity in the
y direction doesn't depend significantly on the cAMP gradient and fluctuates around 4 μm min
^{−1}. 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
^{−1} ≤ ∇
C ≤ 1 nM μm
^{−1} over two orders of magnitude, the chemotactic velocity is almost constant. For gradients above ∇
C = 1 nM μm
^{−1} the directionality of movement is decreased. Exceeding an upper threshold of ∇
C_{up} ∼ 10
^{2} nM μm
^{−1}, the cell motion becomes isotropic again.

 Fig. 3 Independent data analysis carried out by M. Theves^{13} showing (left) chemotaxis as a function of gradient steepness ∇C: above a threshold at ∇C ≃ 10^{−3} nM μm^{−1} cells show a positive (in our coordinate system negative) average velocity in the gradient direction (v_{y}) as well as an increased total motility v, while the perpendicular velocity component in the flow direction (v_{x}) remains random and averages to zero within error bars. For gradients ranging over two orders of magnitude, 10^{−2} nM μm^{−1} ≤ ∇C ≤ 1 nM μm^{−1}, the chemotactic speed is constant. (right) Average chemotactic velocity v_{y} 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 an increase in the average midpoint concentration. This effect seems to reverse for steep gradients above 1 nM μm^{−1}, where v_{y} decreases slightly in higher concentration backgrounds. The figures are used by courtesy of M. Theves from his master's thesis.^{13}  
• 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^{−1}, v_{y} increases as the background concentration increases. This effect reverses for steep gradients ∇C > 0.3 nM μm^{−1} (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 the 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 the 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 previously, chemotactic movement of the cells depends 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 a Langevintype 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, the midpoint concentration is globally averaged out. Other quantities such as chemotactic index, defined as the distance moved in the gradient direction divided by the total distance, are also averaged quantities where information on midpoint concentration is lost. However in our analysis, instead of velocity and acceleration, we work directly with the 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 checking 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} and _{0}, respectively, to obtain 
〈y〉(t) = 〈y〉_{0} + ṽ_{0}t,  (15) 

σ_{y}^{2}(t) = σ_{y}^{2}(0) + 2_{0}t.  (16) 
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 σ_{y}^{2} (see Fig. 2 and 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 over time. This is a motivation to use drift and diffusion parameters which do depend on the position. There is an obvious positiondependence in our system: C (the concentration of cAMP) depends on the coordinate y. Assuming that the drift and diffusion parameters do depend on C, one is left with a ydependence in the drift and diffusion parameters. A simple manageable model is to Taylorexpand this ydependence 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, with wildtype 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 shows a superdiffusive behavior. This superdiffusive 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 datadriven model with purely diffusive behavior. As we mentioned above, a pure diffusive model cannot explain the nonlinear behavior observed in the experimental data of ref. 7.
In our analysis, we observed that at all concentrations D_{0} + D_{1}〈y〉 decrease with time and become 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^{−1} (C_{min} = 0, C_{max} = 50 nM). Comparison between Fig. 2 and 4 shows a similar decreasing behavior in σ_{y}^{2}(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 eqn (21), σ_{y}^{2}(t) = σ_{y}^{2}(0) exp(2v_{1}t). If everything is expanded up to first order in v_{1}, then the result is σ_{y}^{2}(t) = σ_{y}^{2}(0)(1 + 2v_{1}t). As v_{1} ≤ 0, it seems that σ_{y}^{2}(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 positiondependence 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. The fact that the time dependence of the variance matches the experiments, means that the method works. But the perturbative parameters should not be misleading.

 Fig. 4 (a) Trajectories of three different experiments recorded at the same cAMP gradient of 0.14 nM μm^{−1} are combined to improve the statistics. (b) 207 (out of 1097) trajectories, are selected and (c) truncated based on the criteria in Section II.D. These trajectories participate in our analysis which is more than two times the number of selected trajectories in Fig. 2. (d and e) The comparison between experimental data (red lines) and the fitted model (blue line) for 〈y〉 and σ_{y}^{2}.  
Furthermore, with our model, we can test directly the spacetime 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 Fig. 2, 6, 7 and 9–11 seem to show a drift in the positive x direction (in addition to the chemotactic drift in the −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 σ_{x}^{2}(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 conditions. This nonpure shift in σ_{x}^{2}(t) could be a hallmark of correlation between the displacement of individual cells along the 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 a flow speed of 320 μm s^{−1} 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 meanfield approximation.

 Fig. 5 (a) The same trajectories as in Fig. 2 which are shifted to x = 0, (b) selected and (c) truncated based on the criteria in Section II.D. The y dependency of all trajectories is kept as before because the midpoint concentration is different along the width of the channel. (d and e) show the behavior of σ_{x}^{2} 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.  

 Fig. 6 In an experiment with C_{max} = 1 nM (∇C = 0.003 nM μm^{−1}), out of 282 trajectories shown in panel (a), after applying the selection conditions of Section II.D, only 41 trajectories are selected in panel (b) and truncated in panel (c). The comparisons between the experimental data (red) and the model (blue) are presented in panels (d and e).  

 Fig. 7 (a–c) Based on our selection criteria in Section II.D, only 27 (out of 815) trajectories participate in our analysis for C_{max} = 10 nM (∇C = 0.028 nM μm^{−1}). (d and e) show 〈y〉 and σ_{y}^{2} plotted for both experimental data (red) and the model (blue).  

 Fig. 8 (a–c) 61 (out of 650) trajectories satisfy the selection conditions of Section II.D for a gradient of ∇C = 0.09 nM μm^{−1} (C_{max} = 31.6 nM). In panels (d and e) the outcome of experimental data and the model are compared.  

 Fig. 9 Out of 2321 trajectories in panel (a), only 25 are selected in panel (b) and truncated in panel (c) for C_{max} = 100 nM (∇C = 0.28 nM μm^{−1}). A large number of trajectories either show immobile cells or become discontinuous as the cells collide. Again panels (d and c) are the comparison between the data and the model.  

 Fig. 10 (a–c) 39 trajectories (out of 254) are participating in our analysis for C_{max} = 316 nM (∇C = 0.9 nM μm^{−1}). y and σ_{y}^{2} are calculated from the truncated trajectories (red lines) and compared with the fitted model (blue lines) in panels (d and e).  

 Fig. 11 (a–c) In an experiment with C_{max} = 10000 nM (∇C = 28.6 nM μm^{−1}), out of 401 trajectories, 41 are selected and truncated. Comparisons between the experimental measured quantities (red) and the fitted model (blue) are shown in panels (d and e).  
To improve our statistics, we have divided long mother trajectories into 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. has shown that (see the right panel of Fig. 3) with an increase in the average midpoint concentration the average chemotactic velocity v_{y} doesn't show any clear trend for intermediate gradients, 10^{−2} nM μm^{−1} ≤ ∇C ≤ 0.3 nM μm^{−1}. However, for shallow gradients v_{y} increases with midpoint concentration and for steep gradients it decreases. Most of our analysis is carried out at intermediate gradients where the variation in chemotactic velocity v_{y} between three different bins is less than 25 percent. For example, at ∇C = 0.3 nM μm^{−1}, the average chemotactic velocity changes from ∼3 μm min^{−1} to ∼4 μm min^{−1} for three bins with different midpoint concentrations of 17 nM, 50 nM and 83 nM. At steep gradients larger than 0.3 nM μm^{−1}, the variations in v_{y} are even less than 10 percent. Thereby we believe that shorter daughter trajectories which belong 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 will 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 a linear gradient of cAMP. In order to have reliable statistics, we kept the number of trajectories during our analysis constant. Trajectories were selected based on two criteria: (i) they should persist for at least 20 min, and (ii) within this time interval, the cells should have migrated more than 20 μm in the direction of 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 superdiffusive behavior can dominate at larger time scales. In particular, while in a conventional advection–diffusion model both the mean and the variance are linear over time, here in both cases terms arise which are quadratic over time.
In future studies, 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 compositions.^{28} As the substrate composition changes, the cells regulate forces generated by the actomyosin network to maintain the optimal cellsurface 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 where 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 wildtype and mutant cells may reflect important information about the role of the various proteins in the signaling pathway of D.d. cells, which possibly cannot be resolved in models where the midpoint concentration information is 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.
Conflicts of interest
There are no conflicts to declare.
Appendix I
Here we show how to derive eqn (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 
 (17) 
We can directly obtain an equation for the time evolution of 〈y(t)^{2}〉 by multiplying the master equation, eqn (7), by y^{2} and integrating over y. This results in 
 (18) 
The lefthand side of this equation is simply equal to . We apply partial integration to the righthand side and obtain 
 (19) 
Since σ_{y}^{2}(t) ≡ 〈y(t)^{2}〉 − 〈y(t)〉^{2}, and , one finds 
 (20) 
where according to eqn (10), d〈y〉/dt has been replaced by v_{0} + D_{1} + v_{1}〈y(t)〉. The solution of eqn (20) is found as 
 (21) 
where σ_{y}^{2}(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 
σ_{y}^{2}(t) = σ_{y}^{2}(0) + 2[σ_{y}^{2}(0)v_{1} + D_{0} + D_{1}〈y〉_{0}]t + (2D_{0}v_{1} + D_{1}v_{0})t^{2}.  (22) 
The mean displacement of chemotactic cells, in both the x and y directions, and the corresponding variances can be calculated from the experimental data as follows

 (23) 

 (24) 

 (25) 

 (26) 
where
N denotes the number of cells.
Appendix II
Here we show how shifting the cells' tracks to x = 0 can affect on the variance of the xcomponent through the time, σ_{x}^{2}(t). Let x_{i}(t) be the position of the i'th particle in the xdirection at time t. The displacement of the i'th particle in the xdirection through time is z_{i}(t) = x_{i}(t) − x_{i}(0). Simply, one has 〈z(t)〉 = 〈x(t)〉 − 〈x(0)〉 where 〈x(t)〉 and 〈x(0)〉 denote the mean values of the xcomponent of the particles at time t and t = 0, respectively, and 〈z(t)〉 is the mean displacement of the particles. It is worth mentioning that for example , where N denotes the number of cells. In order to find the variance, first we note that 
x_{i}(t) − 〈x(t)〉 = [z_{i}(t) − 〈z(t)〉] + [x_{i}(0) − 〈x(0)〉].  (27) 
After squaring both sides of the above equation and averaging, one finds 
〈[x(t) − 〈x(t)〉]^{2}〉 = 〈[z(t) − 〈z(t)〉]^{2}〉 + 〈[x(0) − 〈x(0)〉]^{2}〉 + 2〈[x(0) − 〈x(0)〉][z(t) − 〈z(t)〉]〉.  (28) 
The covariance of z(t) and x(0) is defined as cov[z(t),x(0)] ≡ 〈(z(t) − 〈z(t)〉)(x(0) − 〈x(0)〉)〉. This quantity provides a measure for the strength of the correlation between two stochastic variables. Using the definition of the variance and covariance, eqn (28) can be written as 
σ_{x}^{2}(t) = σ_{z}^{2}(t) + σ_{x}^{2}(0) + 2cov[z(t),x(0)].  (29) 
We see that when z(t) and x(0) are independent, one has 〈z(t)x(0)〉 = 〈z(t)〉〈x(0)〉 and cov[z(t),x(0)] becomes zero. In this case, σ_{x}^{2}(t) differs from σ_{z}^{2}(t) by just a constant shift. In other words, the necessary and sufficient condition for a pure shift in the variances of σ_{x}^{2}(t) and σ_{z}^{2}(t) is vanishing of the covariance of z(t) and x(0).
Appendix III
As we discussed in the main text, the experiments had been carried out at different cAMP concentrations. Here we present the trajectories, before and after a selection procedure, as well as the corresponding analysis for the cAMP concentrations of C_{max} = 1,10, 31.6, 100, 316, and 10000 nM.
Acknowledgements
We are deeply grateful to M. Theves, E. Bodenschatz, G. Amselem and C. Beta for sharing the experimental data of ref. 7 and A. Bae for critical reading of the manuscript. Z. E. and F. M.R. are grateful to A. Celani, R. Golestanian and L. MollazadehBeidokhti for helpful discussions. Z. E. and F. M.R. acknowledge the hospitality of ICTP in Trieste, where some parts of this work were carried out. F. M.R. acknowledges the hospitality of MPIDSLFPN Group in Göttingen, where this work was initiated. M. K. acknowledges the support of the research council of Alzahra University. A. G. acknowledges the support of the MaxSynBio Consortium which is jointly funded by the Federal Ministry of Education and Research of Germany and the Max Planck Society. Open Access funding provided by the Max Planck Society.
References
 H. Takagi, M. J. Sato, T. Yanagida and M. Ueda, PLoS One, 2008, 3, e2648 Search PubMed .
 W. F. Loomis, D. Fuller, E. Gutierrez, A. Groisman and W. J. Rappel, PLoS One, 2012, 7, 42033 Search PubMed .
 M. Buenemann, H. Levine, W. J. Rappel and L. M. Sander, Biophys. J., 2010, 99, 50 CrossRef CAS PubMed .
 K. F. Swaney, C. H. Huang and P. N. Devreotes, Annu. Rev. Biophys., 2010, 39, 265 CrossRef CAS PubMed .

R. H. Kessin, Dictyostelium: evolution, cell biology, and the development of multicellularity, Cambridge University Press, 2001 Search PubMed .
 C. S. Patlak, Bull. Math. Biophys., 1953, 15, 311 CrossRef CAS .
 G. Amselem, M. Theves, A. Bae, E. Bodenschatz and C. Beta, PLoS One, 2012, 7, e37213 CAS .
 G. Amselem, M. Theves, A. Bae, C. Beta and E. Bodenschatz, Phys. Rev. Lett., 2012, 109, 108103 CrossRef PubMed .

A. Friedman and C.Y. Kao, Mathematical Modeling of Biological Processes, Springer, 2014 Search PubMed .
 E. F. Keller and L. A. Segel, J. Theor. Biol., 1970, 26, 399 CrossRef CAS PubMed .
 H. Othmer and A. Stevens, SIAM J. Appl. Math., 1997, 57, 1044 CrossRef .
 P. Chavanis, Eur. Phys. J. B, 2008, 62, 179 CrossRef CAS .

M. Theves, Quantitative Study of Eukaryotic Chemotaxis with Microfluidic Devices, Master's thesis, GeorgAugust Universität Göttingen, and Max Planck Institut für Dynamik und Selbstoganisation, 2009 Search PubMed .
 N. Andrew and R. H. Insall, Nat. Cell Biol., 2007, 9, 193 CrossRef CAS PubMed .
 L. Bosgraaf and P. J. M. Van Haastert, PLoS One, 2009, 4, e6842 Search PubMed .
 P. J. M. Van Haastert, Biophys. J., 2010, 99, 3345 CrossRef CAS PubMed .
 L. Song, S. M. Nadkarni, H. U. Bodeker, C. Beta and A. Bae,
et al.
, Eur. J. Cell Biol., 2006, 85, 981 CrossRef CAS PubMed .
 L. Li, E. C. Cox and H. Flyvbjerg, Phys. Biol., 2011, 8, 046006 CrossRef PubMed .
 N. Makarava, S. Menz, M. Theves, W. Huisinga, C. Beta and M. Holschneider, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 042703 CrossRef PubMed .
 N. L. Jeon, S. K. W. Dertinger, D. T. Chiu, I. S. Choi, A. D. Stroock and G. M. Whitesides, Langmuir, 2000, 16, 168311 CrossRef .
 E. Decave, D. Rieu, J. Dalous, S. Fache, Y. Brechet, B. Fourcade, M. Satre and F. Bruckert, J. Cell Sci., 2003, 116, 4331 CrossRef CAS PubMed .
 J. C. Crocker and D. G. Grier, J. Colloid Interface Sci., 1996, 179, 298 CrossRef CAS .

J. D. Murray, Mathematical Biology, Springer, Berlin, 1991 Search PubMed .
 D. Selmeczi, S. Mosler, P. H. Hagedorn, N. B. Larsen and H. Flyvbjerg, Biophys. J., 2005, 89, 912 CrossRef CAS PubMed .
 M. Khorrami, A. Shariati, A. Aghamohammadi and A. H. Fatollahi, Phys. Lett. A, 2012, 376, 687 CrossRef CAS .
 P. Dieterich, R. Klages, R. Preuss and A. Schwab, Proc. Natl. Acad. Sci. U. S. A., 2008, 105, 459 CrossRef CAS PubMed .
 A. Upadhyaya, J.P. Rieu, J. A. Glazier and Y. Sawada, Physica A, 2001, 293, 549 CrossRef .
 C. P. McCann, E. C. Rericha, C. Wang, W. Losert and C. A. Parent, PLoS One, 2014, 9, 87981 Search PubMed .

This journal is © The Royal Society of Chemistry 2017 