Modelling of Dictyostelium discoideum movement in a linear gradient of chemoattractant

Abstract

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 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 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–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.⁴ When chemotactic attraction prevails over diffusion, the chemotaxis can trigger a self-accelerating process until aggregation takes place. As a result, 10⁵–10⁶ 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 long-range dispersal.⁵

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.⁶ 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⁶ and Keller¹⁰ 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.¹¹ This leads to ordinary mean field Fokker–Planck equations for cell density with space and time dependent coefficients.¹² 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.¹⁸ 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.¹⁵ The measured probabilities of pseudopod directions were used to obtain an analytical model for chemotaxis of cell populations.¹⁶ The prediction of the model is similar to the 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 and 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 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⁻¹ 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 KH₂PO₄, 1 mM Na₂HPO₄) at a density of 2 × 10⁶ 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 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 of v̄ = 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⁻³ 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

∇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 non-linear area were excluded from statistics. Moreover, given the dimension of the channel and the dynamic viscosity of the flowing phosphate buffer (η = 10⁻³ Pa s), one can calculate the shear stress applied on the cells at the imposed mean flow velocity of v̄ = 320 μm s⁻¹ to be σ = 0.038 Pa. According to the literature, mechanosensing in D.d. cells has been observed above a threshold of σ = 0.5 Pa.²¹ 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.¹³

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²² 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 time-points 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⁻¹). 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².

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 v̄_y > 1 μm min⁻¹. Cells with v̄_y < 1 μm min⁻¹ 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 v̄_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 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.²³ 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 x̂ and ŷ, where x̂ 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 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

J = −D∂_yP + vP. (2)

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

∂_tP = −∂_yJ, (3)

where ∂_t denotes differentiation with respect to t. Using eqn (2) and (3), one can find the diffusion–advection equation for the problem as

∂_tP = ∂_y(D∂_yP) − ∂_y(vP). (4)

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

v = v₀ + v₁y + … (5)

and

D = D₀ + D₁y + … (6)

We keep the terms only up to the first order of y, and treat v₁, and D₁ as perturbation coefficients. Here, we assume that the current in the y direction does not depend on x. Using eqn (4)–(6), one finds

∂_tP = (D₀ + D₁y)∂_y²P + (D₁ − v₀ − v₁y)∂_yP − v₁P. (7)

The mean value of the y-component of the cells' positions is obtained by

Differentiating the above expression with respect to time results

After substituting ∂_tP(y,t) from eqn (7) and integrating, one can find

By solving this simple ordinary differential equation we findwhere 〈y〉₀ ≡ 〈y〉|_t=0 denotes the mean initial y-position of the cells. As it has been mentioned above, v₁ and D₁ 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₁ and D₁, one can findIt is worth mentioning that since terms like v₁D₁ are the second order of perturbation parameters, these terms are dropped.

The variance of the cells' positions along y is defined as σ_y²(t) ≡ 〈y(t)²〉 − 〈y(t)〉². Using a similar method (see Appendix I for details), one can find σ_y²(t) as

σ_y²(t) = σ_y²(0) + 2[σ_y²(0)v₁ + D₀ + D₁〈y〉₀]t + (2D₀v₁ + D₁v₀)t², (13)

where σ_y²(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₀, v₁, D₀ and D₁. 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 σ²_y,exp. Tables 1–3 include the best fit values of the above mentioned parameters at different cAMP gradients. In Table 1, 〈y〉₀ and σ_y²(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	10 880	5454.2
10	0.03	347.22	25 893	6857.1
31.6	0.09	417.88	24 136	7608.9
50	0.14	365.99	27 253	6707.3
100	0.29	410.24	32 453	4719.2
316	0.9	371.13	30 854	7033.1
10 000	28.6	382.33	27 665	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 + v1〈y〉0	v 0 + v1〈y〉0 + D1	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
10 000	−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 + D1〈y〉0	2D0v1 + D1v0	2(σy(0)^2v1 + D0 + D1〈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
10 000	−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₀v₁/2 in eqn (12) is small for all concentrations (see the last column of Table 2).

• The mean square displacement function σ_y²(t) shows decreasing behavior at C_max = 10, 50, 316 and 10 000 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 mid-point 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²(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₀ + D₁〈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 parameters 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₀ and v₁ but also on D₁. The coefficient v₀v₁ 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⁻¹. This is consistent with an independent data analysis performed by M. Theves¹³ (see Fig. 3): within a plateau, ranging from 10⁻² 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² nM μm⁻¹, the cell motion becomes isotropic again.

Fig. 3 Independent data analysis carried out by M. Theves¹³ showing (left) chemotaxis as a function of gradient steepness ∇C: above a threshold at ∇C ≃ 10⁻³ nM μm⁻¹ 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⁻² nM μm⁻¹ ≤ ∇C ≤ 1 nM μm⁻¹, 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⁻¹, where v_y decreases slightly in higher concentration backgrounds. The figures are used by courtesy of M. Theves from his master's thesis.¹³

• For all gradients, while the cells crawl up the gradient, the magnitude of 〈v_y〉 = v₀ + v₁〈y〉 decreases as the midpoint concentration increases. However, the independent data analysis of M. Theves¹³ shows a transition: for shallow gradients, right after the onset of chemotaxis ∇C = 0.003 nM μm⁻¹, v_y increases as the background concentration increases. 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 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.,⁷ 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 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, 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 mid-point 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 ṽ₀ and D̃₀, respectively, to obtain

〈y〉(t) = 〈y〉₀ + ṽ₀t, (15)

σ_y²(t) = σ_y²(0) + 2D̃₀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² (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 position-dependence 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 y-dependence in the drift and diffusion parameters. A simple manageable model is to Taylor-expand 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, with wild-type and mutated epithelial canine kidney cells, it has been shown that the cell dynamics can be characterized by an anomalous diffusion.²⁶ In particular, mean squared displacement shows a super-diffusive behavior. This super-diffusive behaviour was also observed in the mean square displacement of Hydra cells.²⁷ 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 cannot explain the non-linear behavior observed in the experimental data of ref. 7.

In our analysis, we observed that at all concentrations D₀ + D₁〈y〉 decrease with time and become negative for concentrations of 10, 31.6, 316, and 10 000 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 4 shows a similar decreasing behavior in σ_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₀ = D₁ = 0. Then, according to eqn (21), σ_y²(t) = σ_y²(0) exp(2v₁t). If everything is expanded up to first order in v₁, then the result is σ_y²(t) = σ_y²(0)(1 + 2v₁t). As v₁ ≤ 0, it seems that σ_y²(0) could become negative after a while. But that is an artifact of the approximation. We intended to find position-dependence 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. 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⁻¹ 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².

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 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²(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 non-pure shift in σ_x²(t) could be a hallmark of correlation between the displacement of individual cells along the x direction and their initial x-positions (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⁻¹ 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.

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 mid-point concentration is different along the width of the channel. (d and e) show the behavior of σ_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.

Fig. 6 In an experiment with C_max = 1 nM (∇C = 0.003 nM μm⁻¹), 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⁻¹). (d and e) show 〈y〉 and σ_y² 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⁻¹ (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⁻¹). 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⁻¹). y and σ_y² 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 = 10 000 nM (∇C = 28.6 nM μm⁻¹), 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⁻² 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 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⁻¹, the 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 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 super-diffusive 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.²⁸ As the substrate composition changes, the cells regulate forces generated by the actomyosin network to maintain the 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 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 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 cannot be resolved in models where the mid-point 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)²〉 is defined asWe can directly obtain an equation for the time evolution of 〈y(t)²〉 by multiplying the master equation, eqn (7), by y² and integrating over y. This results inThe left-hand side of this equation is simply equal to . We apply partial integration to the right-hand side and obtainSince σ_y²(t) ≡ 〈y(t)²〉 − 〈y(t)〉², and , one findswhere according to eqn (10), d〈y〉/dt has been replaced by v₀ + D₁ + v₁〈y(t)〉. The solution of eqn (20) is found aswhere σ_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₁ and D₁, one has

σ_y²(t) = σ_y²(0) + 2[σ_y²(0)v₁ + D₀ + D₁〈y〉₀]t + (2D₀v₁ + D₁v₀)t². (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

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 x-component through the time, σ_x²(t). Let x_i(t) be the position of the i'th particle in the x-direction at time t. The displacement of the i'th particle in the x-direction 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 x-component 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)〉]²〉 = 〈[z(t) − 〈z(t)〉]²〉 + 〈[x(0) − 〈x(0)〉]²〉 + 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²(t) = σ_z²(t) + σ_x²(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²(t) differs from σ_z²(t) by just a constant shift. In other words, the necessary and sufficient condition for a pure shift in the variances of σ_x²(t) and σ_z²(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 10 000 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. Mollazadeh-Beidokhti 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 MPIDS-LFPN 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

1. H. Takagi, M. J. Sato, T. Yanagida and M. Ueda, PLoS One, 2008, 3, e2648.
2. W. F. Loomis, D. Fuller, E. Gutierrez, A. Groisman and W. J. Rappel, PLoS One, 2012, 7, 42033.
3. M. Buenemann, H. Levine, W. J. Rappel and L. M. Sander, Biophys. J., 2010, 99, 50.
4. K. F. Swaney, C. H. Huang and P. N. Devreotes, Annu. Rev. Biophys., 2010, 39, 265.
5. R. H. Kessin, Dictyostelium: evolution, cell biology, and the development of multicellularity, Cambridge University Press, 2001.
6. C. S. Patlak, Bull. Math. Biophys., 1953, 15, 311.
7. G. Amselem, M. Theves, A. Bae, E. Bodenschatz and C. Beta, PLoS One, 2012, 7, e37213.
8. G. Amselem, M. Theves, A. Bae, C. Beta and E. Bodenschatz, Phys. Rev. Lett., 2012, 109, 108103.
9. A. Friedman and C.-Y. Kao, Mathematical Modeling of Biological Processes, Springer, 2014.
10. E. F. Keller and L. A. Segel, J. Theor. Biol., 1970, 26, 399.
11. H. Othmer and A. Stevens, SIAM J. Appl. Math., 1997, 57, 1044.
12. P. Chavanis, Eur. Phys. J. B, 2008, 62, 179.
13. M. Theves, Quantitative Study of Eukaryotic Chemotaxis with Microfluidic Devices, Master's thesis, Georg-August Universität Göttingen, and Max Planck Institut für Dynamik und Selbstoganisation, 2009.
14. N. Andrew and R. H. Insall, Nat. Cell Biol., 2007, 9, 193.
15. L. Bosgraaf and P. J. M. Van Haastert, PLoS One, 2009, 4, e6842.
16. P. J. M. Van Haastert, Biophys. J., 2010, 99, 3345.
17. L. Song, S. M. Nadkarni, H. U. Bodeker, C. Beta and A. Bae, et al. , Eur. J. Cell Biol., 2006, 85, 981.
18. L. Li, E. C. Cox and H. Flyvbjerg, Phys. Biol., 2011, 8, 046006.
19. N. Makarava, S. Menz, M. Theves, W. Huisinga, C. Beta and M. Holschneider, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 042703.
20. N. L. Jeon, S. K. W. Dertinger, D. T. Chiu, I. S. Choi, A. D. Stroock and G. M. Whitesides, Langmuir, 2000, 16, 168311.
21. E. Decave, D. Rieu, J. Dalous, S. Fache, Y. Brechet, B. Fourcade, M. Satre and F. Bruckert, J. Cell Sci., 2003, 116, 4331.
22. J. C. Crocker and D. G. Grier, J. Colloid Interface Sci., 1996, 179, 298.
23. J. D. Murray, Mathematical Biology, Springer, Berlin, 1991.
24. D. Selmeczi, S. Mosler, P. H. Hagedorn, N. B. Larsen and H. Flyvbjerg, Biophys. J., 2005, 89, 912.
25. M. Khorrami, A. Shariati, A. Aghamohammadi and A. H. Fatollahi, Phys. Lett. A, 2012, 376, 687.
26. P. Dieterich, R. Klages, R. Preuss and A. Schwab, Proc. Natl. Acad. Sci. U. S. A., 2008, 105, 459.
27. A. Upadhyaya, J.-P. Rieu, J. A. Glazier and Y. Sawada, Physica A, 2001, 293, 549.
28. C. P. McCann, E. C. Rericha, C. Wang, W. Losert and C. A. Parent, PLoS One, 2014, 9, 87981.

---