Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Phys Rev E Stat Nonlin Soft Matter Phys. Author manuscript; available in PMC 2010 May 19.
Published in final edited form as:
Phys Rev E Stat Nonlin Soft Matter Phys. 2010 March; 81(3 Pt 1): 031906.
Published online 2010 March 8. doi:  10.1103/PhysRevE.81.031906
PMCID: PMC2872938

Phenomenological approach to eukaryotic chemotactic efficiency


Eukaryotic cells are capable of detecting small chemical gradients for a wide range of background concentrations. Ultimately, fluctuations place a limit on gradient sensing and recent work has focused on the role of stochastic receptor occupancy as one possible limiting factor. Here, we use a phenomenological approach to add spontaneous motility fluctuations to receptor noise and predict the directional statistics of eukaryotic chemotaxis. Specifically, an Itô diffusion equation with direction-dependent multiplicative noise is developed and analytically studied. We show that our approach can naturally accommodate recent experimental data for the chemotaxis of the social amoeba Dictyostelium.


During chemotaxis, cells direct their movement by sensing chemical gradients. Chemotaxis plays an important role in various biological processes including fertilization, neuronal development, wound healing and cancer metastasis [1, 2]. Unlike prokaryotes, eukaryotic cells use a spatial measurement of the asymmetric distribution of bound receptors following the external gradient. These bound receptors trigger the activation of intracellular signaling pathways, eventually leading to the generation of directional movement. Interestingly, eukaryotic cells are able to direct their motion in shallow gradients. For example, the social amoeba Dictyostelium discoideum is observed to crawl up a chemoattractant gradient when the front-back difference in the concentration is around 1% [3, 4]. Furthermore, the cells respond even when the average local concentration is well below the dissociation constant Kd. These findings has led to significant interest from the biophysics community in understanding the effect of noise on chemotaxis [512].

Studies to date have mainly focused on fluctuations in gradient sensing, arising for example from stochastic receptor dynamics [6, 10]. However, there are other sources of noise, as is evidenced by the fact that cells move randomly in the absence of directional information [13, 14]. The simplest approach to adding random motility to directional bias is to assume that the cell executes Brownian motion in a deterministic periodic potential [15]; however this approach fails to take into account uncertainties in evaluating gradient steepness and direction. Thus, what is needed is a new model of directed cell motility that takes into account both sources of fluctuations. Such a model should be able to relate the relative importance of gradient-determination noise versus motility noise to testable predictions for the behavior of chemotactic cells. Here, we propose a phenomenological, stochastic differential equation (SDE) model for this purpose and compare its predictions with the available data. As we will see, our approach naturally accounts for the unusual shape of the directional probability distribution for chemotaxing cells.


In our model, the cellular surface takes a circular geometry where N independent receptors are uniformly distributed (Fig. 1). We divide the surface into M sectors such that the local ligand concentration near each sector is constant while the number of receptors in each sector Nm = N/M is large enough for the continuous approximation. For Dictyostelium, a typical value is N = 60000, and one may choose M = 300 and Nm = 200. We define the gradient steepness parameter as p=LC0dCdr, where C0 is the average local ligand concentration across the cell length L. Then, the local concentration at the mth sector with angular position [var phi]m is given by Cm=C0[1+p2cos(φmϕ)] where ϕ specifies the gradient direction. The time for receptor decor-relation is dominated by the individual receptor rate τm = (k + k+Cm)−1 which is typically faster than the process by which motion can be mechanically altered [16]. This allows us to use the white-noise approximation in which the number of occupied receptors is Ym = left angle bracketYmright angle bracket + ηm = NmCm/(Cm + Kd) + ηm for m = 1, …, M, with left angle bracketηm(tn(s)right angle bracketNmCmKd/(Cm + Kd)2 δ(tsmn [7, 17]. Thus, the receptor signal is decomposed into M independent Gaussian random variables.

(Color online). Schematic representation of our model. The forward rate k+ and backward rate k determine the transition between the unoccupied R0 and occupied R1 states. The dissociation constant Kd is the ligand concentration at which half the ...

The simplest estimate a cell can make of the gradient is obtained by the following statistic: Z=m=1MYmcosφm+im=1MYmsinφmz1+iz2. When M is large, one can evaluate Z via replacing the sum by an integral. For small gradients (p < 10%), the integrand can be expanded around p and we find



It is easy to check that z1 and z2 are uncorrelated Gaussian random variables, with different means but the same variance. This implies that Z is a complex Gaussian variable which can be written in polar coordinates as Z = ρeiψ. Here, ρ measures the degree of asymmetry in the ligand-bound receptor distribution while the phase ψ estimates the gradient direction ϕ. Mathematically, ρ follows the Rice distribution and ψ takes a symmetric unimodal circular distribution [18]. It was found that for large signal-to-noise ratio (SNR= ν/σ [greater, similar] 3) both ρ and ψ are asymptotically Gaussian [18]. Thus, ρ ≈ ν + ηψ with left angle bracketηψ(tρ(s)right angle bracket = σ2δ(ts), and ψ ≈ ϕ + ηψ with left angle bracketηψ(tψ(s)right angle bracket = (σ22)δ(ts). As an orthogonal transformation from the Cartesian coordinates, ηρ is independent of ηψ. For typical Dictyostelium values we find that SNR> 3, consistent with the Gaussian approximation.

Cellular motion can be decomposed into two stochastic processes: speed and direction. Recent experimental data in Dictyostelium [4] show that the average speed is roughly identical for all directions. In these experiments, cells were exposed to stable exponential chemoattractant profiles using microfluidic devices. The centroid of each cell in the field of view was automatically tracked every 5s for 100 frames. Cells that moved the furthest without colliding with another cell were chosen for data analysis. Ten to twenty-five such cells were collected in each experiment, giving thousands of data points for each particular gradient. The average speed as a function of the angle is plotted in Fig. 2 for two different gradient steepnesses. From this figure, we can conclude that any possible correlation between angle and speed is smaller than can be obtained from the data. Thus, we restrict ourselves to the directional process, parametrized by the migration angle θ(t), by assuming a uniform migration speed.

The average cell speed in the experiments of Ref. [4] as a function of angle for p=5% (A) and p=10% (B). For each gradient steepness, the speed was computed by tracking the position of the centroid of 30–40 cells during 8 minutes. For more experimental ...

Let us temporarily suppose that cells have perfect knowledge about the gradient direction ϕ. In general, the equation of θ(t) can be written as dθ/dt = G(ρ, θ − ϕ) + η0 along with η0(t)η0(s)=σ02δ(ts). The incorporation of the gradient-independent noise η0 above allows the cell to perform a random walk in the absence of any gradient. Note that unlike a recent study of random motility [13], we do not include colored noise in our model. The function G(ρ, θ) must exhibit the symmetry properties [15]: (1) 2π-periodicity, G(ρ, θ ± 2nπ) = G(ρ, θ) for any integer n; (2) polar symmetry, G(ρ, θ ± π) = −G(ρ, θ), because the system’s polarity will be reversed under the reversal of gradient direction; (3) reflection symmetry, G(ρ, −θ) = −G(ρ, θ), i.e. cells cannot distinguish between left and right with respect to the direction ϕ. The simplest form obeying these requirements is G(ρ, θ − ϕ) = −β(ρ) sin (θ − ϕ). Of course, cells have only imperfect knowledge of the gradient direction. The phase variable ψ in our model serves as an unbiased estimator of ϕ for the cell. Thus, we arrive at the phenomenological Langevin equation


This phenomenological equation describes a chemotactic cell trying to align its movement with the estimated gradient direction ψ. The first term describes a gradient-dependent restoring force which steers the cell toward the gradient direction and which contains gradient-dependent fluctuations (arising from receptor-ligand dynamics and downstream pathways) while the second term represents gradient independent noise in the random motility machinery.

Without loss of generality, we will set ϕ = 0. In the small noise limit we can perform a Taylor expansion dθ/dt ≈ −β(ν) sin θ + β′(ν) sin θηρ + β(ν) cos θηψ + η0 = −β(ν) sin θ + ηtot, where the variance of the total noise is


with β′(ν) = [partial differential]β(x)/[partial differential]x|x and β−1(ν) = β(ν)/ν. The first term in Eq. (4) describes the fluctuations arising from ρ, the amplitude of the gradient estimate, while the second term describes the fluctuations in ψ, the direction of the gradient estimate. Both of them vary with the instantaneous direction of motion θ(t) and therefore are multiplicative fluctuations. Since the cell responds to noise at its current position, the multiplicative noise is defined with the Itô prescription. From the previous expression it is clear that the ultimate directional dependence in σtot2 relies on the sign of β′(ν) − β−1(ν). To capture possible nonlinear input-output characteristics of downstream pathways, we assume that β(ν) has a power-law functional form, i.e. β(ν) = aνb. We can distinguish between two qualitatively different cases. In the first, b > 1, the total variance is of the “sin θ” type, σtot2=σ02+(σβ1)2+(β2β12)(σsinθ)2 and is maximal at θ = π/2; in the second, b < 1, the variance is described by a “cos θ” type and is maximal at θ = 0. In the special case b = 1, the total variance is reduced to σtot2=σ02+a2σ2 with no directional dependence.


We first study the case of b > 1 where the restoring force is ultra-sensitive to the input signal. Define the direction-independent noise σ2=σ02+(σβ1)2 and the direction-dependent noise σ2=(β2β12)σ2. Then the previous Langevin equation can be approximated as the Itô SDE


where η(t) is white noise. One can solve the associated Fokker-Planck equation on the interval [−π, π] by imposing periodic boundary conditions [19]. The stationary distribution reads


where κ=2β/σ2 is a parameter that compares the direction-biased force with the direction-independent noise and where λ=σ/σ2+σ2 is the ratio of the direction dependent and maximal total variance. The normalization parameter Ω can be determined numerically. For a vanishing λ, we recover the Circular Normal (CN) distribution: limλ→0 Ps(θ; κ, λ) = exp(κ cos θ)/(2πI0(κ)). For large λ, our stationary distribution will appreciably deviate from the CN distribution with a sharper peak at θ = 0 and heavier tails near θ = ±π (see Fig. 3A). A physical interpretation of this leptokurtic feature can be given by realizing that the first term in the Langevin equation (Eq. (3)) functions as a restoring force. Contrary to the CN case, the amplitude of this force, β(ρ), is taken from a distribution. Thus, large values can occur and correspond to large restoring forces which lead to a sharper peak in Ps(θ). Using the same argument, one can see that the occurrence of small values in this distribution are responsible for the heavy tails in Ps(θ).

The stationary distribution Ps(θ|ϕ = 0; κ = 1.5, λ) for (A) b > 1 and (C) b < 1 as a function of the parameter λ representing the ratio of direction-dependent noise to maximal total variance. The ...

We can further explore the extrema and directionality of Ps(θ) by recalling that arctanh (x)=12In(1+x1x) and setting εκ(1λ2)/(2λ)=λβ(ν)/σ2, after which we obtain the alternative form Ps(θ) = Ω(1+λ cos θ)ε−1/(1 − λ cos θ)ε+1. The solution of [partial differential]θ ln Ps = 0 contains five possible roots in the interval [−π, π]: θ = 0, θ = ±π, and θ = ± arccos(−ε/λ) which exist only if λ > ε. As seen in Fig. 3A, Ps(θ = ±π) correspond to global minima if λ < ε, but become local maxima once λ > ε, in which case the new global minima are located at θ = ± arccos(−ε/λ). This change of extrema occurs at a critical point λc = ε (equivalent to σ2=β) where the direction-dependent noise becomes comparable to the magnitude of the restoring force. Direct numerical simulations of the full model described by Eq. (3) indicate that our analytical expression for Ps(θ) works well when SNR> 3 and λ < ε, justifying the Gaussian noise approximation.

As is common in the experimental literature, we can define an order parameter left angle bracketcos θright angle bracket, the Chemotaxis Index (CI), which is a measure of the directionality of the angular distribution. In Fig. 3B, we plot the CI for different values of κ using our model. The dashed line in Fig. 3B corresponds to the CI with a value of κ that satisfies the critical condition λ = ε or equivalently κ = 2λ2/(1 − λ2). To the right of this line the distribution shows a local maximum at ±π. The slope of the CI is close to 0 over a large range of values for λ. This can be understood by realizing that the negative effect of the heavy tails on the CI is partially compensated by the positive effect of a sharpened peak.

Next we look at the case of b < 1. The corresponding stationary distribution is


with σ2=σ02+(σβ)2 and σ2=(β12β2)σ2 in the definitions for κ and λ. The analysis proceeds as above and we only highlight the differences. The angle distribution exhibits not only heavy tails but also an obtuse peak for increasing values of λ (Fig. 3C). This can be understood since in this case the noise in ψ, i.e. the second term in Eq. (4), dominates and reduces the accuracy in directionality. Due to the Gaussian approximation, the analytical probability distribution exhibits double peaks at θ = ± arccos(−ε/λ) when λ > ε [20]. Also, the CI decays quickly with increasing λ even before the bifurcation (Fig. 3D), suggesting that cells should operate in a regime where b ≥ 1. Finally, for b = 1, the directional dependence in σtot2 disappears and the stationary distribution reduces to the CN distribution with shape parameter κ˜=2aν/(σ02+a2σ2). Obviously, for larger values of the parameter a, the system is more sensitive to the gradient, and more of the receptor noise is transmitted to the directional output.

For any given a and σ02, there exists a unique value of b that maximizes the CI (Fig. 4A). We found numerically that the combination of a* = σ0/σ and b* = 1 optimizes the CI over the whole parameter space (a, b) given the level of σ0. The optimal pair (a*, b*) corresponds to the ridge of the CI surface in Fig 4A. In summary, the exponent b is critical for the chemotactic responses. The optimal output requires b = 1, while b > 1 is suboptimal and b < 1 leads to poor directionality. This can be understood by realizing that when restoring force β(ν) is a nonlinear function (i.e. b ≠ 1) of the receptor signal ν, the associated sensor noise is modulated to be direction-dependent, which leads to heavy tails in the angle distribution (Fig. 3). In other words, there is a significant probability that cells may move in the wrong direction and the chemotactic performance is not optimal. An important implication of our model is that the CI is always lower than that predicted by models neglecting receptor noise (see Fig. 4B). In fact, previous experimental data indicates that when the gradient steepness increases the saturating value of the CI is less than one [5, 9]. This feature is not consistent with previous CN models but may be successfully explained by our b > 1 model where the total noise σtot2 becomes an increasing function of p.

(Color online). (A) CI maximized by choosing appropriate b for different values of a and σ02, given p = 5% and C0 = Kd = 30nM. (B) CI vs. p for C0 = Kd = 30nM, σ02=2, b = 1, and a = a*, 2.5a/*, and 0.25a*. The solid line corresponds to ...

We can compare our model to the results of recent experiments in which Dictyostelium cells were exposed to stable exponential chemoattractant profiles [4]. The data allow us to compute directional distributions for different gradient parameters. To eliminate the slight asymmetric bias introduced by the flow in the microfluidic chambers we collected the data into 20 bins according to their absolute deviation from the gradient direction, |θ−ϕ|. The resulting distribution are shown as symbols for p = 5% and p = 10% in Fig. 4C and 4D, respectively. We have also plotted as a solid line a least squares fit to the data using our model (Eq. (6)). Our fit captures well the heavy tails and the sharp peaks, features that the CN distribution, plotted as a dashed line, is unable to reproduce. A goodness of fit analysis using a Pearson’s chi-square test revealed that the CN distribution is rejected at a significance level of less than 1% while our model exhibits large p-values [21]. The observed leptokurtic feature suggests that Dictyostelium cells are operating close to the optimal regime with the exponent b estimated as slightly larger than 1.


A key assumption in our model is that the chemotaxing cell experiences a relatively stable gradient steepness (p = const) regardless of the cell’s location or moving history. Thus, our model results can be directly compared to experiments in which the gradient is carefully quantified and controlled. If the gradient steepness is spatially varying, the chemotactic response would be heterogeneous in space, possibly obscuring the structure of the angular distribution.

Although we focus here on the role of receptor noise, it is easy to extend our model to include extra gradient-dependent fluctuations, by adding them to the magnitude of the restoring force, β(ρ). Such fluctuations might arise from noise in the downstream signal transduction pathways. Future work will track cells for longer periods of time and test for the presence of long-lived cell individuality and other fluctuations sources [22]. Should the data indicate such memory effects, our model could be extended to include quenched fluctuations in cell responses.


DF, WFL, HL and WJR were supported by NIH Grant P01GM078586. BH acknowledges support from NSF PHY-0822283. We thank Y.-H. Tu, R. J. Williams, J. Wolf, and T. Hwa for useful discussions.


PACS numbers: 02.50.Le, 05.65.+b, 87.23.Ge, 87.23.Kg


1. Haastert PJV, Devreotes PN. Nat. Rev. Mol. Cell Biol. 2004;5:626. [PubMed]
2. Parent CA, Devreotes PN. Science. 1999;284:765. [PubMed]
3. Song L, Nadkarni SM, Bödeker HU, Beta C, Bae A, Franck C, Rappel W-J, Loomis WF, Bodenschatz E. Eur. J. Cell Biol. 2006;85:981. [PubMed]
4. Fuller D, et al. (submitted)
5. Haastert PJV, Postma M. Biophys. J. 2007;93:1787. [PubMed]
6. Wang K, Rappel W-J, Kerr R, Levine H. Phys. Rev. E. 2007;75:061905. [PubMed]
7. Rappel W-J, Levine H. Phys. Rev. Lett. 2008;100:228101. [PMC free article] [PubMed]
8. Rappel W-J, Levine H. Proc. Natl. Acad. Sci. U.S.A. 2008;105:19270. [PubMed]
9. Endres RG, Wingreen NS. Proc. Natl. Acad. Sci. U.S.A. 2008;105:15749. [PubMed]
10. Bialek W, Setayeshgar S. Proc. Natl. Acad. Sci. U.S.A. 2005;102:10040. [PubMed]
11. Shibata T, Fujimoto K. Proc. Natl. Acad. Sci. U.S.A. 2005;102:331. [PubMed]
12. Ueda M, Shibata T. Biophys. J. 2005;93:11. [PubMed]
13. Li L, Nørrelykke SF, Cox EC. PLoS ONE. 2008;3:e2093. [PMC free article] [PubMed]
14. Takagi H, Sato MJ, Yanagida T, Ueda M. PLoS ONE. 2008;3:e2648. [PMC free article] [PubMed]
15. Schienbein M, Franke K, Gruler H. Phys. Rev. E. 1994;49:5462. [PubMed]
16. Here we have used the fact that under typical conditions the contribution of ligand diffusion to the receptor correlation can be ignored [6, 7].
17. Lauffenburger DA, Liderman JJ. Receptors: Models for Binding, Trafficking, and Signaling. New York: Oxford University Press; 1993.
18. Gudbjartsson H, Patz S. Magn. Reson. Med. 1995;34:910. [PMC free article] [PubMed]
19. Risken H. Fokker-Planck Equation. Berlin: Springer-Verlag; 1984.
20. Andrews BW, Iglesias PA. PLoS Comput. Biol. 2007;3:e153. [PubMed]
21. For p = 5%, the chi-square test of our model results in χour2=14.14, with a degree of freedom DF = 20 − 3 = 17 and a p-value equal to 0.657. Here, the DF is reduced by 3 since the model has two fitting parameters and the number of data points is constrained. For comparison, the CN distribution gives χCN2=81.54 with DF = 18 and p-value< 10−3. For p = 10%, we find χour2=10.51, p-value= 0.881, and χCN2=694.78, p-value< 10−3.
22. Samadani A, Mettetal J, van Oudenaarden A. Proc. Natl. Acad. Sci. U.S.A. 2006;103:11549. [PubMed]