On the efficiency of chemotactic pursuit - Comparing blind search with temporal and spatial gradient sensing

In chemotaxis, cells are modulating their migration patterns in response to concentration gradients of a guiding substance. Immune cells are believed to use such chemotactic sensing for remotely detecting and homing in on pathogens. Considering that immune cells may encounter a multitude of targets with vastly different migration properties, ranging from immobile to highly mobile, it is not clear which strategies of chemotactic pursuit are simultaneously efficient and versatile. We tackle this problem theoretically and define a tunable response function that maps temporal or spatial concentration gradients to migration behavior. The seven free parameters of this response function are optimized numerically with the objective of maximizing search efficiency against a wide spectrum of target cell properties. Finally, we reverse-engineer the best-performing parameter sets to uncover strategies of chemotactic pursuit that are efficient under different biologically realistic boundary conditions. Although strategies based on the temporal or spatial sensing of chemotactic gradients are significantly more efficient than unguided migration, such ‘blind search’ turns out to work surprisingly well, in particular if the immune cells are fast and directionally persistent. The resulting simulated data can be used for the design of chemotaxis experiments and for the development of algorithms that automatically detect and quantify goal oriented behavior in measured immune cell trajectories.


For all cells c >=0 (including the immune cell c=0):
Draw step width w c,n from Rayleigh distribution, according to momentary speed parameter v c,n of cell c.
Draw turning angle ΔΦ c,n from uniform distribution, according to momentary persistence parameter ε c,n of cell c.
Compute magnitude |ΔΦ c,n | of turning angle.
Draw sign s c,n of turning angle from Bernoulli distribution, according to momentary 'right turn probability' q R;c,n of cell c (Note: q R;c,n = ½ for all target cells c>0).
Compute new cell position r c,n , according to Eq.(2).

For all cells c >= 0 (including the immune cell c=0) :
Apply periodic boundary conditions.

For all target cells c>0:
Compute distance d c0 to immune cell 0, taking into account periodic boundary conditions. Update previous value ρ C n-1 at the central chemoattractant sensor. Update present value ρ C n at the central chemoattractant sensor. Compute temporal gradient Δρ C n , according to Eq.(15). Compute momentary 'approach mode probability' q A;0,n , according to Eq.(18).

Do with probability q A;0,n :
Toggle immune cell speed v 0,n between v N and v A . Toggle immune cell persistence ε 0,n between ε N and ε A .
Update value ρ R n at the right chemoattractant sensor. Update value ρ L n at the left chemoattractant sensor. Compute spatial gradient Δρ LR n , according to Eq.(17). Compute momentary 'right turn probability' q R;0,n , according to Eq.(20). In this paper, we focused on the special case of a chemoattractant that is diffusing fast (large D) and decaying quickly (large k), so that the critical velocity v crit = √ kD becomes much larger than the typical migration speed of the target cells. In this so-called 'fast diffusion limit', the global density profile F 2D ( r, t) could be computed by the linear superposition of fixed 'kernels' f 2D ( r), one centered around each of the target cells.
In the general case, it is necessary to simulate the temporal evolution of the density field, which is governed by Eq. (6), along with the motion of the immune and target cells.
The main structure of the simulation algorithm (Supplemental Information [A]) remains the same, but the subroutine 'Update density field' becomes more elaborate: The most simple way to solve the partial differential equation Eq. (6) is based on a discretization of the simulation area into a regular grid of quadratic patches with linear patch size ∆L, so that the continuous density field F 2D ( r, t) is represented by the discrete array F 2D (X, Y, n). In every simulation time step n, the chemoattractant density in each of the patches is changing due to the generation term, the diffusion term, and the decay term in Eq. (??). The generation term effectively adds, in every time step, a certain amountĝ of density to all patches (X, Y ) in which a target cell is presently located. The diffusion term removes an amountD F 2D (X, Y, n) from every patch (X, Y ), and redistributes this amount in equal parts to the four nearest neighbor patches (X+1, Y ), (X, Y +1), (X−1, Y ) and (X, Y −1). Finally, the decay term multiplies the density at every patch (X, Y ) with a certain factorm that is smaller than one. Given the spatial discretization length ∆L and the temporal discretization time ∆t sim , the three effective simulation parametersĝ,D and m can be computed from the actual model parameters g, D and k. However, we will use the effective parameters in the following two examples.
In the examples, we consider target cells that move diffusively ( tar = 0) and at low speed (v tar = 1). The immune cell moves with the same speed (v imm = 1), but uses temporal gradient sensing to switch its directional persistence between imm = 0.7 in the approach mode and imm = 0 in the normal migration mode. Response coefficients were c A0 = 10 and c A1 = 500. The effective parameter for chemoattractant generation wasĝ = 1, and the effective diffusion was (D = 0.02).
If the chemoattractant decays relatively quickly (m = 0.9, Fig. 3), the situation is similar to the fast diffusion limit: Each target cell is surrounded by a distinct 'cloud', and the gradient within these clouds guides the immune cell reliably to the targets themselves.
A qualitatively different situation arises if the chemoattractant decays relatively slowly (m = 0.999, Fig. 4). In this case, the emissions of several targets 'bleed together' and form flat density plateaus. It is difficult for the immune cells to pick up density gradients within these plateaus.
These example demonstrate that the efficiency Q of different search strategies (such as TGS and SGS) will probably depend in a complex way on the diffusion constant D and decay rate k of the chemoattractant. While it would be interesting to map out the complete 'phase diagram' Q = Q(v tar , tar , D, K), this is clearly beyond the scope of this paper.