Self-Generated Chemoattractant Gradients: Attractant Depletion Extends the Range and Robustness of Chemotaxis

Chemotaxis is fundamentally important, but the sources of gradients in vivo are rarely well understood. Here, we analyse self-generated chemotaxis, in which cells respond to gradients they have made themselves by breaking down globally available attractants, using both computational simulations and experiments. We show that chemoattractant degradation creates steep local gradients. This leads to surprising results, in particular the existence of a leading population of cells that moves highly directionally, while cells behind this group are undirected. This leading cell population is denser than those following, especially at high attractant concentrations. The local gradient moves with the leading cells as they interact with their surroundings, giving directed movement that is unusually robust and can operate over long distances. Even when gradients are applied from external sources, attractant breakdown greatly changes cells' responses and increases robustness. We also consider alternative mechanisms for directional decision-making and show that they do not predict the features of population migration we observe experimentally. Our findings provide useful diagnostics to allow identification of self-generated gradients and suggest that self-generated chemotaxis is unexpectedly universal in biology and medicine.


Computational model
The model consists of a diffusible attractant and a population of cells at discrete locations, able to degrade attractant locally. The profile of attractant concentration c(x, t) changes as follows: r(c(x n, t , t))δ(x, x n, t ) (1) Here the first term represents the diffusion of the attractant, with a diffusion constant D c . The second term sums the contributions of a population of N cells to the degradation of attractant, with the nth cell found at the location x n, t at time t. Each cell degrades attractant explicitly locally (enforced here by the Kronicker delta) and at a rate r, which depends on the local concentration. r takes the form: where r max is the rate of decay when the cell is saturated with attractant, and k M is the Michaelis-Menten constant.
Practically, diffusion of attractant in the model was simulated on a grid using a central differences approximation of the diffusion equation. The degradation term was simulated using the Euler method. Cells do not move on the diffusion grid, and so can take intermediary positions. As such, degraded attractant was taken from surrounding grid points, linearly interpolating the rate based on their distances from the cell. In simulations of the Insall chamber, the boundary conditions c(0, y, t) = 0 and c(L, y, t) = c max are applied, approximating the effects of the large resevoirs of medium at either end of the bridge.
Cells in simulations of the one-well assay begin in uniformly-distributed random positions in a thin (500µm) band along the left-hand side of the simulation. They move from their positions at a constant speed in a persistent, biased random walk. In each time step of size ∆t, the nth cell chooses a direction θ n, t in which it will move in the interval (t, t + ∆t], by taking a weighted circular average of a persistent random walk p n, t and a bias introduced by local chemical cues b n, t (c(x, t), x n, t ): The direction of the persistent, random contribution to the cell's movement is a random variable chosen from a wrapped normal distribution, centred on the previous direction of motion θ n, t−∆t .
where σ controls the strength of the persistence, and can be seen as the standard deviation in direction per unit time. For the bias b n, t (c(x, t), x n, t ) we choose the difference in the probability of occupancy for a receptor between the very front and back of the cell. Assuming the cell is projected across some length l: The change in the position of the cell over the period (t, t + ∆t] is then ∆x n, t = v c cos(θ n, t )î + sin(θ n, t )ĵ ∆t.
Contact inhibition of locomotion is defined as occuring when two cells come into contact, repolarise and move in a new direction. We simulate this by giving each cell a region of occupancy, and testing other cells for contact when they move. If a contact is detected, the component of movement in the direction of the contacted cell is removed. That is, for cell m, which has made contact with cell n, the new update in position ∆x * m is 2 Mathematical models of the population density wave

Population density wave in a semi-unbounded domain
Let us consider a concentration distribution of attractant c(x, t), initially evenly distributed at a concentration c 0 for all x > 0 which evolves according to the diffusion equation: This attractant is degraded rapidly by a population of cells, which act as a perfect sink. These cells begin in a well at x = 0. As the well is long and the cells evenly distributed, we can treat this as a one dimensional problem. Let the population of cells move over time to a position s(t), which is reached by simple chemotaxis, advancing proportionally to the local gradient: Eqn. (8) is satisfied by the function where erf x denotes the error function of x. As c( where erfc x is the complementary error function of x. Substituting this into Eqn (10) gives us As the wave at s(t) is a perfect sink, we know that c(s(t), t) = 0. Thus, We can use Eqns (17, 18) to estimate the value of α for our system from the behaviour of the wave. Fitting the mean x position of forward-moving cells over time, we estimate that α = 1.23. We can use Eqns. (9) and (12) to show that: and so can use our estimate of α to calculate a value for k. Though the wave moves over time, attractant only reaches the wave over minuscule length scales by diffusion. We can calculate the flux of attractant across the wave using Fick's first law: The minus sign here only serves to tell us that the flux comes from a higher x position, and so we can ignore it. How many cells must the wave contain in order to degrade all the attractant it encounters? If each cell can degrade attractant at a maximum rate r, and the the number of cells required is N we can equate rN to the flux at the position of the wave,

Decay of attractant profile in a bounded domain
Here we consider the profile of a chemoattractant, c(x, t), stretching from a well containing cells at x = 0 to the end of the domain at x = L. We will assume that cells move very slowly relative to other processes, such that we need not update their position. The profile of the chemoattractant again changes according to Eqn (8), the diffusion equation.
We use variable separation in order to arrive at the Fourier series solution for this equation. We assume that the cell population is able to degrade all attractant that reaches it by diffusion very rapidly. This gives us the boundary condition c(0, t) = 0, restricting us to only the Fourier sine series. We also impose a zero-flux boundary condition ∂c(L,t) ∂x = 0, restricting us to half-integer terms. Each remaining term in the series decays exponentially in time. The solution is: where α n = (n + 1/2)π/L. To determine the values of our A n , we impose the initial condition of a square wave of amplitude c 0 . This gives us a uniform profile, save for at x = 0.
The total flux of chemoattractant onto the cells is again given by Fick's first law Or, equivalently, this can be written using the Jacobi theta function As with the unbounded case, we can then relate the flux and the required population N with showing us again that the minimum number of cells needed in order to degrade all attractant decays over time. Here the real position for each frame is taken to be the position of cells moving away from the well, weighted by the speed of this directed movement. The model value for α s fitted to these data. (R) Attractant profile predicted by the simple model after 5h. The gradient shows us why we can expect to see a density increase in the wave, even without saturation-guidance cues are best at the back of the wave, where all attractant is degraded. Ahead of this point, guidance cues are poorer due to the lower gradient, and so the wave will catch up with any cells that get ahead.

Comparison of analytical and experimental results
In these analytical models, the parameter α controls the movement of the wave, such that when α = 0, the wave is stationary at x = 0. When α → 0 and L is large there is (unsurprisingly) good agreement between the two models. Indeed, for α = 0, L = 30mm and D c = 12, 000µm 2 /s, the required number of cells predicted by each model only differs by 0.25% over a simulation time of 5 days (though at this stage both models are, of course, unrealistic). As the contribution of α to the population of the wave happens entirely through the constant 1 erfc α e −α 2 , over moderate time-scales the effect of migration is simply to increase the required population by a constant factor. Fig. S1 shows the population, wave position and attractant profile predicted by the semi-unbounded model. Wave position is calculated by the average position of cells moving away from the well, weighted by their outward velocity. The model represents the population and position reasonably well, though the position of the wave does not match well at short times. We attribute this to limits on the speed and persistence of real cells not present in this idealised analytical model. The predicted gradient also demonstrates why the population density wave forms; cells ahead of the wave position experience a shallower gradient and by extension worse guidance cues. This means that they progess more slowly than the back of the wave does, and fall back into it.

The chemotactic fraction of a domain
A system has a point source of chemoattractant at x = 0 that maintains its concentration at c 0 . It also has a sink at x = L, such that the steady state solution for the system is The general solution here is the sum of two solutions. One is the steady state solution, Eq. (11), the other is a Fourier series. For the latter, we apply the boundary conditions c(0, t) = c(L, t) = 0, restricting us to the integer harmonics of a Fourier sine series. We then find the corresponding Fourier components by equating it with the initial condition of the steady state solution minus a sawtooth wave of wavelength 2L, The general solution we obtain is The expected receptor occupancy κ(x, t) for a receptor at position x and t is and so the relative occupancy difference between the extreme ends of a cell of length l, centred at position x > l/2 is: This is of course a generous simplification, as most receptors will not be at these extrema. We can use Eqs. 12 and 14 to calculate ∆κ(x, t) for any x and t. Zigmond gives the threshold value for chemotaxis as a relative difference in occupancy of 1%. We use this threshold to determine the fraction of the domain which can provide chemotactic guidance, the fraction of x ∈ (0, L) for which ∆κ(x, t) ≥ 0.01. As c(x, t) and k d only appear in the ratio c(x, t)/k d it is clear that what determines guidance is the source concentration relative to receptor affinity, not either of these quntities alone. Similarly, as D c and t only appear in the product D c t, faster diffusion and longer waiting times are equivalent for transporting guidance cues from a point source. As such, we use these expressions as our axes for Fig. 8. Our numerical calculations assume that l = 10µm.