Modelling of Dictyostelium discoideum movement in linear gradient of chemoattractant

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~\mu m/min$ on glass substrates Takagi-2008 ; Loomis:2012 ; Sander-2010 . The chemotactic cell motion is highly organized over a length scale significantly larger than the size of a single cell ($\sim$10 $\mu$m). When nutrients are depleted, D.d. cells secret a chemical called cAMP (cyclic adenosine monophosphate) that has an attractive effect on the cells themselves. Cells sense gradients of cAMP and direct their chemotactic movements towards regions of higher concentration of cAMP Devreotes-2004 . When chemotactic attraction prevails over diffusion, the chemotaxis can trigger a self-accelerating process until aggregation takes place. As a result, $10^{5}$ – $10^{6}$ cells stream towards the aggregation centers and eventually transform into millimeter long slugs and ultimately form fruiting bodies bearing spores for long-term survival and long-range dispersal book-Kessin .

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 Patlak-1953 . 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 book-Friedman-Kao ; Segel-1970 ; Amselem-2012 ; Amselem-2012-PRL . 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 Patlak-1953 and Keller Segel-1970 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 Othmer-1997 . This leads to ordinary mean field Fokker-Planck equations for cell density with space and time dependent coefficients Chavanis-2008 . On the other hand, if we assume that the diffusion coefficient and the drift velocity of the cells depend on their concentration and on the concentration of the secreted chemical, the original Keller-Segel model takes the form of a generalized mean field Fokker-Planck equation.

The statistical characteristics of trajectories of motile D.d. cells have been the subject of several recent studies. These include experiments to characterize chemotactic cell movement in homogenous and inhomogenous chemical cues Amselem-2012 ; Amselem-2012-PRL ; Theves:2009 ; Song-2006 ; Andrew-2007 ; Haastert-2009 ; Haastert-2010 ; Flyvbjerg-2011 , and parallel theoretical modeling to reproduce statistical features of the experimental observations Amselem-2012 ; Holschneider:2014 ; Song-2006 ; Haastert-2010 ; Flyvbjerg-2011 . Recently, Li et al. have presented an experimental study of the individual cells in a homogeneous medium Flyvbjerg-2011 . 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 Haastert-2009 . The measured probabilities of pseudopod directions were used to obtain an analytical model for chemotaxis of cell populations Haastert-2010 . The prediction of the model are similar to measured chemotactic index of wild-type cells as well as the mutants. Besides, although it is well-known that the directed movement of the D.d. cells in response to the chemoattractant cAMP depends both on the absolute value of the local concentration (chemokinesis) and its gradient (chemotaxis), the exact dependency is not well understood.

In this study, we aim to extract the concentration dependencies of the diffusion and drift coefficients in the Fokker-Planck equation (with respect to cell density), by analyzing the experimental trajectories of motile D.d. cells in Ref. Amselem-2012 ; Amselem-2012-PRL ; Theves:2009 . We assume that these coefficients depend on both the local cAMP concentration (the so-called midpoint concentration) and its gradient. The experiments are performed in a microfluidic device (see Section II-B) that generate linear stable gradients between the two inlet concentrations $C_{\text{max}}$ and $C_{\text{min}}$. As the cells crawl up the gradient, the average background concentration they experience increases. These experiments systematically explore different steepnesses and cover a wide range of gradients, at which chemotactic behavior is observed. In these experiments, an external flow removes the naturally produced cAMP secreted by the cells to avoid cell-cell signaling. This is completely different from aggregation process where the cell density is much higher and the cells signal each other. We start our analysis by assuming linear dependencies for the diffusion coefficient and the drift velocity of the cells along the width of microfluidic setup, where a linear gradient is established. We then use the experimental cell trajectories to deduce the coefficients of these linear dependencies at different cAMP gradients.

Nonlinear mean field Fokker-Planck equations can find important applications in the context of chemotaxis Murray-1991 . Here, we attempt to implement an advection-diffusion approach to describe the chemotactic movement of the D.d. cells experiencing a linear stationary gradient Amselem-2012 ; Amselem-2012-PRL . 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 $\hat{x}$ and $\hat{y}$, where $\hat{x}$ is the flow direction and $-\hat{y}$ is the direction of the spatial gradient of cAMP (see Fig. 2). The position of each cell is given by $\vec{r}=x\hat{x}+y\hat{y}$. The concentration of the D.d cells is low enough, so we can assume that each cell does not sense the presence of the other cells. As stated in Experiment section, microfluidic gradient along $y$ direction is generated by a continuous flow along $x$. To avoid mixing up the issues of chemotaxis in response to the chemoattractant and mechanotaxis under the influence of the shear stress due to viscous forces, we limit our model to the chemotactic movement of the cells along $y$. Let us assume that $p(x,y,t)$ denotes the number density of cells at position $(x,y)$ at time $t$. Then, we have the probability density $P(y,t)$, which is the original density $p(x,y,t)$ integrated over $x$. The current density along $y$, $J$, reads as

where $v$ and $D$ are drift velocity and diffusion coefficient, respectively, and $\partial_{y}$ means differentiation with respect to $y$. Now, the continuity equation for $P$ and $J$ then reads as

where $\partial_{t}$ denotes differentiation with respect to $t$. Using Eqs. (2) and (3), one can find the diffusion-advection equation for the problem as

The chemotactic motion of the cells depends on both the absolute local concentration (chemo-kinesis) and its gradients (chemotaxis) Amselem-2012 ; Song-2006 . Based on the experiments of Ref. Amselem-2012 , here we consider a constant spatial gradient of cAMP in the direction of $-\hat{y}$. Therefore, one can expect that both the diffusion coefficient, $D$, and the drift velocity, $v$, depend on the $y$ component of the position vector. Since there is no direct experimental method to determine this dependency, it is plausible to expand the mentioned coefficients in terms of $y$ as

and

We keep the terms only up to the first order of $y$, and treat $v_{1}$, and $D_{1}$ as perturbation coefficients. Here, we assume that the current in the $y$ direction does not depend on $x$. Using equations (4) to (6), one finds

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

Differentiating the above expression with respect to time results

After substituting $\partial_{t}P(y,t)$ from Eq. (7) and integrating, one can find

By solving this simple ordinary differential equation we find

where $\langle y\rangle_{0}\equiv\langle y\rangle|_{t=0}$ denotes the mean initial $y$-position of the cells. As it has been mentioned above, $v_{1}$ and $D_{1}$ are the small parameters and in our model, they have been considered as perturbation parameters. After expanding the exponential factor and keeping the terms up to the first order of perturbation parameters, $v_{1}$ and $D_{1}$, one can find

It is worth mentioning that since terms like $v_{1}D_{1}$ are the second order of perturbation parameters, these terms are dropped.

The variance of the cells’ positions along $y$ is defined as $\sigma_{y}^{2}(t)\equiv\langle y(t)^{2}\rangle-\langle y(t)\rangle^{2}$. Using similar method (see Appendix I for details), one can find $\sigma_{y}^{2}(t)$ as

where $\sigma_{y}^{2}(0)$ is the initial variance of the cells’ positions along $y$. We note that in Eq. (13), we have kept the terms up to the first order of perturbation parameters as well.

Now we are in a position to determine the perturbation parameters of our model based on the experimental trajectories. The mean displacement of chemotactic cells and the corresponding variance can be calculated from the experimental trajectories as defined in Appendix I (Eqs. 23–26). To characterize the chemotactic behavior of D.d. cells, based on Eqs. (12) and (13), we need to determine the values of $v_{0}$, $v_{1}$, $D_{0}$ and $D_{1}$. We treat these factors as tuning parameters and find their values simultaneously by fitting (MATLAB, R2016b) the model to the experimental values of $\langle y\rangle_{\text{exp}}$ and $\sigma_{y,\text{exp}}^{2}$. Tables 1-3 include the best fit values of the above mentioned parameters at different cAMP gradients. In Table 1, $\langle y\rangle_{0}$ and $\sigma^{2}_{y}(0)$ denote the fitted mean and standard deviations at time zero.

Fig. 2 and Figs. 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 are 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-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, $\langle y\rangle$, decreases almost linearly in 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 Eq. 12 is small for all concentrations (see the last column of Table. 2).

  

• •

  The mean square displacement function $\sigma^{2}_{y}(t)$ shows decreasing behavior at $C_{max}=10$, $50$, $316$ and $10000~nM$. This trend is an experimental observation, independent of the introduced model, and has to do with the fact that the cells tend to migrate to areas with high cAMP concentration (top part of the channel). Since in these areas mid-point cAMP concentration is high and most of the cAMP receptors are saturated, therefore the cells slow down and accumulate at the top of the channel. This “accumulation” can give rise to a decreasing $\sigma^{2}_{y}(t)$. In the other word, the possible decrease in the variance is due to a drift towards the top areas of the channel.

• •

  The diffusion coefficient in $y$ direction $D_{0}+D_{1}\langle y\rangle$ is initially positive for all concentrations but as the cells migrate upwards and the value of $\langle y\rangle$ 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

  $$v_{\text{drift}}=\frac{d\langle y\rangle}{dt}=v_{0}+v_{1}\langle y\rangle_{0}+D_{1}+v_{0}v_{1}t,$$

  (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_{\text{drift}}$ is essentially constant in 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 $y$ direction doesn’t depend significantly on the cAMP gradient and fluctuates around 4 $\mu m$/min. This is consistent with an independent data analysis performed by M. Theves Theves:2009 (see Fig. 3): within a plateau, ranging from $10^{-2}~nM/\mu m\leq\nabla C\leq 1~nM/\mu m$ over two orders of magnitude, the chemotactic velocity is almost constant. For gradients above $\nabla C=1~nM/\mu m$ the directionality of movement is decreased. Exceeding an upper threshold of $\nabla C_{up}\sim 10^{2}~nM/\mu m$, the cell motion becomes isotropic again.

• •

  For all gradients, while the cells crawl up the gradient, the magnitude of $\langle v_{y}\rangle=v_{0}+v_{1}\langle y\rangle$ decreases as the midpoint concentration increases. However, the independent data analysis of M. Theves Theves:2009 shows a transition: for shallow gradients, right after the onset of chemotaxis $\nabla C=0.003~nM/\mu m$, $v_{y}$ increases as the background concentration rises. This effect reverses for steep gradients $\nabla C>0.3~nM/\mu m$ (see right panel of Fig. 3).

We have analyzed large data sets of D.d. chemotaxis in linear gradients of cAMP recorded by Theves et al. in a microfluidic setup Amselem-2012 ; Amselem-2012-PRL ; Theves:2009 . Data sets with different steepnesses of cAMP gradient were included in our analysis, covering a large range of gradients, in which chemotactic behavior was observed. Inspired by the experimental conditions of Ref. Amselem-2012 , we introduced a minimal phenomenological model that explicitly incorporated the dependency of diffusion matrix and velocity of the cells on their positions which corresponds to the position dependence of the local concentration of chemotactic cues. Based on this model, we extracted the physical properties of the chemotactic D.d. cells using the mean and variance of the experimental cell tracks. What is the benefit of this phenomenological model? As highlighted previously, chemotactic movement of the cells depend on both the chemoattractant gradient and the average ambient chemoattractant concentration (midpoint concentration). In the microfluidic setup of Theves et al. Amselem-2012 , the cells are exposed to a constant gradient, while the midpoint concentration increases when the cells are moving up the gradient. Traditionally, chemotactic cell motion is described by Langevin-type equation where for each cell track, the velocity and acceleration of the cells are calculated at each point by finite differences from the cell positions Amselem-2012 ; Amselem-2012-PRL ; Flyvbjerg-2005 . Therefore in these types of analysis midpoint concentration is globally averaged out. Other quantities such as chemotactic index, defined as the distance moved in gradient direction divided by the total distance, are also averaged quantities where information on mid-point concentration is lost. However in our analysis, instead of velocity and acceleration, we work directly with spatial position of the cells and explicitly include the dependence of the diffusion coefficient and the drift velocity on the midpoint concentration. Taylor expansion of these coefficients up to the first order in $y$ leads to a closed set of equations that can be solved to obtain the fitting parameters. It is worth to check the effect of the dependency of the diffusion and velocity on the local concentration. To do this, let us assume that the drift velocity and the diffusion coefficient were constant. We denote the constant drift velocity and diffusion constant by $\tilde{v}_{0}$ and $\tilde{D}_{0}$, respectively, to obtain

These equations predict a linear dependency on $t$, in both the mean and the variance of the position, which is not consistent with the experimental data especially in the variance of $\sigma_{y}^{2}$ (see Fig. 2 and Figs. 6–11). In fact it has been shown in Ref. Khorrami-2012 , that any linear diffusion model (even anomalous), which enjoys both time translational invariance and space translational invariance leads to means and variances which are at most linear in time. This is a motivation to use drift and diffusion parameters which do depend on the position. There is an obvious position-dependence in our system: C (the concentration of cAMP) does depend on the coordinate $y$. Assuming that the drift and diffusion parameters do depend on C, one is left with a $y$-dependence in the drift and diffusion parameters. A simple manageable model is to 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. Amselem-2012 fulfill the necessary conditions of the mentioned study.

In previous studies, wild-type and mutated epithelial canine kidney cells, it has been shown that the cell dynamics can be characterized by an anomalous diffusion Dieterich-2008 . 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 Upadhyaya-2001 . However, experimental trajectories of chemotactic D.d. cells in Ref. Flyvbjerg-2011 , were interpreted by a data-driven model with purely diffusive behavior. As we mentioned above, a pure diffusive model can not explain the non-linear behavior observed in the experimental data of Ref. Amselem-2012 .

In our analysis, we observed that at all concentrations $D_{0}+D_{1}\langle y\rangle$ decreases with time and becomes negative for concentrations of $10,~31.6,~316,$ and $10000~nM$. In order to make sure that negative diffusion coefficients are not due to our low statistics after the selection procedure, we combined trajectories of three different experiments performed at the same cAMP gradient, namely $\nabla C=0.14~nM/\mu m~(C_{\text{min}}=0,~C_{\text{max}}=50~nM)$. Comparison between Fig. 2 and Fig. 4 shows a similar decreasing behavior in $\sigma^{2}_{y}(t)$. Indeed, an absorbing point on top of the channel can produce a decreasing variance, not through the diffusion but through the upstream drift. Let us suppose that $D_{0}=D_{1}=0$. Then, according to Eq. 21, $\sigma^{2}_{y}(t)=\sigma^{2}_{y}(0)\exp(2v_{1}t)$. If everything is expanded up to first order in $v_{1}$, then the result is $\sigma^{2}_{y}(t)=\sigma^{2}_{y}(0)(1+2v_{1}t)$. As $v_{1}\leq 0$, it seems that $\sigma_{y}^{2}(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$. That the time dependence of the variance does match the experiments, means that the method works. But the perturbative parameters should not be misleading.

Furthermore, with our model, we can test directly the space-time symmetries of the cell movement. Based on the reports of the experiments the gradient in the $y$ direction is homogeneous in $x$. However, the cell tracks shown in the panel (a) of Figs. 2, 6, 7, 9, 10, and 11 seem to show a drift in positive $x$ direction (in addition to the chemotactic drift in -$y$ direction). To check the spatial homogeneity in the $x$ direction, we shifted all the tracks of Fig. 2 to $x=0$ (see Fig. 5). It is interesting that in this case $\sigma_{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 condition. This non-pure shift in $\sigma^{2}_{x}(t)$ could be a hallmark of correlation between the displacement of individual cells along $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 flow speed of $320~\mu m/s$ we are far below the regime, where mechanotactic effects have been observed in D.d. cells. However, the authors of Ref. Mechanotaxis conducted their experiments with vegetative cells. This suggests that starvation may increase the mechanosensitivity of D.d. cells. We emphasize that with this correction all of our analysis in the $y$ direction is still valid, given that the current in the $y$ direction does not depend on $x$. This assumption is nothing but a mean-field approximation.

To improve our statistics, we have divided long mother trajectories to shorter ones and if the criteria in Section. II.D are satisfied, we have included daughter trajectories as completely independent tracks in our analysis. The main difference between these new daughter trajectories is the average midpoint concentration that the cells experience as they crawl up the gradient. This corresponds to moving up from Bin 3 to Bin 1 in Fig. 1d, where in each Bin cells are exposed to a different average midpoint concentration. Detailed analysis by Theves et al. have shown that (see the right panel of Fig. 3) with a raise in the average midpoint concentration the average chemotactic velocity $v_{y}$ doesn’t show any clear trend for intermediate gradients, $10^{-2}nM/\mu m\leq\nabla C\leq 0.3~nM/\mu m$. However, for shallow gradients $v_{y}$ increases with midpoint concentration and for steep gradients it decreases. Most of our analysis are done at intermediate gradients where the variation in chemotactic velocity $v_{y}$ between three different Bins is less than 25 percent. Exemplary, at $\nabla C=0.3~nM/\mu m$, average chemotactic velocity changes from $\sim 3~\mu m/min$ to $\sim 4~\mu 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/\mu m$, the variations in $v_{y}$ is 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 can be helpful to obtain longer trajectories. Experiments with lower cell density (to avoid cell-cell collision) can also help us to obtain longer trajectories, as the cells after collision are indistinguishable from each other and two new trajectories are detected by the cell tracking algorithms.

In summary, we have analyzed the experimental data of chemotactic D.d. cells in the linear gradient of cAMP. In order to have a reliable statistics, we kept the number of trajectories during our analysis constant. Trajectories were selected based on two criteria: (i) they should persist at least 20 min, and (ii) within this time interval, the cells should have migrated more than 20 $\mu$m in the direction of the the gradient of cAMP. We have shown that by introducing an advection-diffusion model that includes the position dependence on the cAMP concentration, a quantitative description of experimental cell tracks of the amoeba D.d. is achieved. Our analysis goes beyond a pure diffusive model and shows that the super-diffusive behavior can dominate at larger time scales. Specifically, while in a conventional advection-diffusion model both the mean and the variance are linear in time, here in both cases terms arise which are quadratic in time.

In future study, we aim to apply our analysis to the trajectories of cells migrating on surfaces of differing composition. In a recent study, it has been shown that D.d. cells migrate similarly on surfaces with various chemical composition Losert:2016 . As the substrate composition changes, the cells regulate forces generated by actomyosin network to maintain optimal cell-surface contact area and adhesion. We will assess migration trajectories of the cells on different surfaces and investigate the variations in the fitting parameters of our model. Furthermore, we aim to extend our analysis to the trajectories of mutant cell lines that single or multiple components of the chemotactic signaling pathway are deficient and consequently the character of the cell trajectories may change considerably. Structural differences between the trajectories of wild-type and mutant cells may reflect important information about the role of the various proteins in the signaling pathway of D.d. cells, which possibly can not be resolved in the models that mid-point concentration informations are averaged out. The objective is to correlate various parameters of our model to the key molecular players involved in chemotaxis. This can provide a link between the observed macroscopic dynamics and the underlying microscopic mechanism which is an important goal in the field of eukaryotic chemotaxis.

We are deeply grateful to M. Theves, E. Bodenschatz, G. Amselem and C. Beta for sharing the experimental data of Ref. Amselem-2012 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 is done. 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 the 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.

Here we show how to derive Eq. (13). Indeed, for what follows we do not need the exact form of $P(y,t)$ itself, but just the time dependence of its moments. $\langle y(t)^{2}\rangle$ is defined as

We can directly obtain an equation for the time evolution of $\langle y(t)^{2}\rangle$ by multiplying the master equation, Eq. (7), by $y^{2}$ and integrate over $y$. This results in

The left-hand side of this equation is simply equal to $\frac{d\langle y(t)^{2}\rangle}{dt}$. We apply partial integration to the right-hand side and obtain

Since $\sigma_{y}^{2}(t)\equiv\langle y(t)^{2}\rangle-\langle y(t)\rangle^{2}$, and $\frac{d}{dt}\sigma_{y}^{2}(t)=\frac{d}{dt}\langle y(t)^{2}\rangle-2\langle y(t)\rangle\frac{d}{dt}\langle y(t)\rangle$, one finds

where according to Eq. (10), $d\langle y\rangle/dt$ has been replaced by $v_{0}+D_{1}+v_{1}\langle y(t)\rangle$. The solution of Eq. (20) is found as

where $\sigma_{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

The mean displacement of chemotactic cells, in both $x$ and $y$ directions, and the corresponding variances can be calculated from the experimental data as follows

where $N$ denotes the number of cells.

Here we show how shifting the cells’ tracks to $x=0$ can affect on the variance of the x-component through the time, $\sigma_{x}^{2}(t)$. Let ${\displaystyle x_{i}(t)}$ be the position of the $i$’th particle in $x$-direction at time $t$. The displacement of the $i$’th particle in $x$-direction through the time is ${\displaystyle z_{i}(t)=x_{i}(t)-x_{i}(0)}$. Simply, one has $\langle z(t)\rangle=\langle x(t)\rangle-\langle x(0)\rangle$ where $\langle x(t)\rangle$ and $\langle x(0)\rangle$ denote the mean values of $x$-component of the particles at time $t$ and $t=0$, respectively, and $\langle z(t)\rangle$ is the mean displacement of the particles. It is worth mentioning that for example $\langle x(t)\rangle\equiv 1/N\sum_{i}x_{i}(t)$, where $N$ denotes the number of cells. In order to find the variance, first we note that

After squaring both sides of the above equation and averaging, one finds

The covariance of $z(t)$ and $x(0)$ is defined as $\mathrm{cov}[z(t),x(0)]\equiv\left\langle(z(t)-\langle z(t)\rangle)(x(0)-\langle x(0)\rangle)\right\rangle$. This quantity provides a measure for the strength of the correlation between two stochastic variables. Using the definition of the variance and covariance, Eq. (LABEL:A3E1) can be written as

We see that when $z(t)$ and $x(0)$ are independent, one has $\langle z(t)x(0)\rangle=\langle z(t)\rangle\,\langle x(0)\rangle$ and $\mathrm{cov}[z(t),x(0)]$ becomes zero. In this case, ${\sigma^{2}_{x}(t)}$ differs from ${\displaystyle\sigma^{2}_{z}(t)}$ by just a constant shift. In other words, the necessary and sufficient condition for a pure shift in the variances of $\sigma_{x}^{2}(t)$ and $\sigma_{z}^{2}(t)$ is vanishing of the covariance of $z(t)$ and $x(0)$.

As we discussed in the main text, the experiments had been done in different cAMP concentrations. Here we present the trajectories, before and after selection procedure, as well as the corresponding analysis for the cAMP concentrations of $C_{\rm max}=1,\;10,~31.6,~100,~316,$ and $10000$ nM.