Principles of efficient chemotactic pursuit
Claus Metzner




claus.metzner@gmail.com
Department of Physics, Biophysics Group, Friedrich-Alexander University Erlangen, Germany




	
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 an 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 takle 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 the principles of efficient chemotactic pursuit under different biologically realistic boundary conditions. Remarkably, the numerical optimization rediscovers chemotactic strategies that are well-known in biological systems, such as the gradient-dependent swimming and tumbling modes of E.coli. Some of our results may also be useful 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.



  Keywords
Chemotaxis, search strategies, immune cells.



Principles of efficient chemotactic pursuit
    Claus Metzner
    December 30, 2023
============================================





§ INTRODUCTION


Chemotaxis, the ability of cells to detect and follow concentration gradients of specific chemicals, is ubiquitus in biology (For an introduction to the field, see <cit.> and the references therein). It helps sperm cells to find the ovum, directs cell movements during embryogenesis, but also enables organisms to locate food sources and to avoid hostile environments. In particular, chemotaxis plays a vital role in recruiting motile immune cells to sites of infection or to malignant tumors. This recruitment of immune cells is often based on endogenous chemo-attractants, which are released by other host cells that are already at the location where a pathogen has invided the body. However, the fact that individual immune cells are able to find and eliminate tumor cells in a Petri dish, without any assistance, suggests that immune cells can also detect chemical traces emitted by the pathogens themselves. 

In this paper, we are mainly interested in this latter scenario of a single immune cell finding and eliminating several target cells on a two-dimensional plane with periodic boundary conditions. The target cells are modeled as simple agents that move, independently from the immune cell and from each other, with fixed speed and with fixed directional persistence. While migrating, the targets are emitting a chemical substance that acts as a chemo-attractant for the immune cell. This chemo-attractant is assumed to spread quickly within the extracellular medium by linear diffusion. It is also assumed to decay at a constant rate, so that a concentration profile of fixed shape will surround each target cell at any moment.  

The immune cell is modeled as a more complex agent with concentration sensors for the chemo-attractant and with the ability to change its migration behavior accordingly. In the simplest case, the immune cell has only a single chemo-attractant sensor and compares the measured local concentrations between subsequent simulation time steps (temporal sensing). In the more powerful case of spatial sensing, the immune cell uses multiple sensors at different body positions to measure the spatial gradient of the chemo-attractant concentration. 

In order to modulate the migration properties depending on the sensed concentration gradients, the immune cell uses  probabilistic 'stimulus-response functions' with tunable parameters. In the case of temporal sensing, the response function controls the momentary probabilities for being in one of two possible modes of migration, characterized by different speeds and degrees of directional persistence. In the case of spatial sensing,
the response function determines the probability of the immune cell turning clockwise or counter-clockwise. 

The parameters of the response functions are optimized numerically, with the objective to maximize the average number of direct contacts between the immune cell and distinct target cells during a fixed simulation time - a number called the 'search efficiency' Q (Here, we assume that once a direct contact is established, the respective target cell is immediately removed from the system). In order to obtain an immune cell that is not only efficient in finding specific types of targets but also robust against variable target behavior, the simulated immune cell is confronted with a broad spectrum of target cell speeds v_tar and directional persistences ϵ_tar during the optimization phase. Once the optimal response parameters are found, we also evaluate the specific performance Q=Q(v_tar,ϵ_tar) of the immune cell as a function of the target cell's migrational properties. 




§ METHODS




 §.§ Cell migration model

We consider a single immune cell (with index c=0) and several target cells (with indices c=1… N_tar) on a two-dimensional simulation area of linear dimension L_sys. The migration of the cells is described by the time-dependent position r⃗_c,n of the respective cell centers, where periodic boundary conditions are applied both in x- and y-direction. Here, n is a discrete time index, related to the continuous time by t_n = n Δ t_sim.


Throughout this work, we use a fixed simulation time interval of

    Δ t_sim := 1   min.



The cell trajectories r⃗_c,n are modelled as discrete time, correlated random walks. In particular, the update from one position to the next is performed as follows:

    r⃗_c,n = r⃗_c,n-1 + w_c,n·(
    [ cos( ϕ_c,n-1+s_c,n |Δϕ_c,n| ); sin( ϕ_c,n-1+s_c,n |Δϕ_c,n| ) ]).



In Eq.(<ref>), w_c,n is the step width, which is randomly and independently drawn from a Rayleigh distribution with mean value v. Note that this corresponds to an average speed of the cell along the contour of the trajectory (which is a sequence of line segments). 


The quantity ϕ_c,n-1 is the planar angle of motion during the last step of cell c, that is, ϕ_c,n-1=arctan(y_c,n-1-y_c,n-2/x_c,n-1-x_c,n-2).


The quantity Δϕ_c,n is the turning angle between the last and the present step of cell c, so that ϕ_c,n = ϕ_c,n-1 + Δϕ_c,n. The turning angles are randomly and independently drawn from a uniform distribution between the limits Δϕ_min(ϵ) and
Δϕ_max(ϵ). Here, ϵ∈[-1,+1] is a persistence parameter, where ϵ=+1 corresponds to fully persistent motion, ϵ=0 to diffusive motion, and ϵ=-1 to fully anti-persistent motion. Consequently, if ϵ>0, we define Δϕ_min(ϵ)=-(1-ϵ)π and Δϕ_max(ϵ)=+(1-ϵ)π. If ϵ<0, we define Δϕ_min(ϵ)=(1-|ϵ|)π and Δϕ_max(ϵ)=(1+|ϵ|)π. Note that only the magnitude of the turning angle enters in Eq.(<ref>)


The quantity s_c,n∈{-1,+1} is a sign factor, which controls if the cell moves left (counter-clockwise) of right (clock-wise). It is randomly and independently assigned to one of its two possible values, with a probability prob(R)=prob(s_c,n = -1)=q_R. 


The statistical properties of the random walk generated by Eq.(<ref>) are determined by the three parameters v, ϵ, and q_R, where v controls the speed of the cells, ϵ their directional persistence, and q_R their preference to turn left or right (which is usually balanced, so that q_R=1/2). In simple cell migration models, these parameters are usually considered as constant over time. However, it has been shown that cell migration is a heterogeneous stochastic process, in which all parameters can change gradually or abruptly, depending on the circumstances of the cell <cit.>. In this work, we assume in particular that the immune cell is able to adapt its speed, persistence and left/right preference in response to local gradients of a chemo-attractant.  



 §.§ Assumed size and migration parameters of cells


If not stated otherwise, simulations in this paper assume that both the immune cell and the target cells are rotation-symmetric and have a radius of

    r_imm=r_tar := 10 μ m.



Target cells are assumed to be slow and to move diffusively:

    v_tar := 1 μ m/min;  ϵ_tar := 0.



If required, the immune cell is able to move much faster than the targets and, at least for short periods, with perfect directional persistence:

    v_imm := 0…6 μ m/min;  ϵ_imm := 0…1.



Indeed, experiments have shown that natural killer cells can migrate in thin collagen gels with an average speed of v_nk≈ 5.67 μ m/min and an average persistence of ϵ_nk≈ 0.64 [Christoph Mark and Franzsika Hoersch, private communication. Note that CM and FH used a slightly different persistent parameter, which however has the same property that ϵ=+1 corresponds to fully persistent motion, ϵ=0 to diffusive motion, and ϵ=-1 to fully anti-persistent motion].



 §.§ Model for temporal evolution of the chemo-attractant


Our basic proposition is that the target cells emit a substance into the extra-cellular matrix (mainly consisting of water), which is used as a chemo-attractant by the immune cell. For simplicity, we assume that the chemo-attractant is produced at the center point r⃗_0 of each target cell with a constant generation rate g. The substance is freely diffusing with diffusion constant D, and is spontaneously decaying with a rate k (It is important  - and also biologically realistic - that this decay rate is non-zero. Otherwise no stationary density profile will develop). This leads to the following partial differential equation for the time-dependent 2D density distribution of the chemo-attractant f_2D(r⃗,t):

    d/dt f_2D = g δ(r⃗-r⃗_0) + D (∇_x^2+∇_y^2)  f_2D - k f_2D,




 §.§ Typical parameters of diffusion and decay


The diffusion constant of a substance within a liquid medium (here basically water) can be estimated by assuming a spherical shape of the diffusing molecules. Using Stokes formula for the friction force, the resulting Stokes-Einstein relation yields

    D = k_BT/6πη r,

where T=37 ^∘C is the temperature, η=6.91·10^-4 Pa s is the  viscosity of water at this temperature, and r is the radius of the diffusing molecule. For a hypothetic molecule with r=3.18  nm, one obtains a diffusion constant of 

    D := 100 μ m^2/s,

which will be used throughout this paper. Note that the same value of D was used in an analytical study of the chemo-attractant's density profile <cit.>, where the considered molecule was the anaphylatoxin C5a. Following this reference, we also assume a typical decay constant of 

    k := 10^-2/s.

The generation rate g is less important in the sense that it does not affect the spatial shape or the temporal evolution of the profile f_2D(r⃗,t).


A dimensional analysis of Eq.(<ref>) reveals that the system has a characteristic diffusion length of 

    L_dif = √(D/k)≈ 100 μ m.
      ,

which can be considered as the approximate spatial extent of the density 'cloud' around a stationary emitter. The characteristic time period for developing this density cloud can be estimated as

    T_dif = L_dif^2/D = 1/k≈ 100 s ≈ 1.7  min.
      ,




 §.§ Fast diffusion limit


Based on the above parameters, we can compute a further characteristic quantity that has the dimensions of a velocity:

    v_crit = L_dif/T_dif≈ 60 μ m/ min
      ,

If the emitter of the density cloud is moving at a speed much smaller than this critical velocity, we can approximately assume that the density cloud is fully developed at any moment in time. In other words, there will be a cloud of fixed (stationary) shape that is 'carried around' by the emitter along its trajectory. For our assumed typical target cell speed of v_tar=1  μ m/s, we are indeed well within this 'fast diffusion limit'.



 §.§ Stationary density profile around single target


The fast diffusion limit saves us from numerically solving the reaction-diffusion equation Eq.(<ref>). We only need to compute the stationary, rotation-symmetric density profile f_2D( r = |r⃗-r⃗_0| ) around a non-moving emitter, conveniently located at the origin r⃗_0=0⃗ of the coordinate system. Since the immune cell can never be closer to the emission point r⃗_0 than the radius r_tar of the target cells, we need to solve Eq.(<ref>) only in the region r>r_tar, where the generation term disappears. In polar coordinates, these simplifications lead to the following ordinary differential equation for the stationary profile,

    [∂^2/∂ r^2 +1/r ∂/∂ r]  f_2D(r) = (k/D)  f_2D(r)
      ,

which is solved numerically with a Runge Kutta method. At the border of the prey cell, without restriction of generality, we set the density to f_2D(r=r_prey)=1. The slope ∂/∂ rf_2D(r=r_prey) at this point is iteratively adjusted such that f_2D(r→∞)=0. The resulting radial profile decays rapidly in the direct vicinity of the emitter. For r→∞, the decay approaches an exponential shape (see Fig.<ref>(b)).



 §.§ Momentary distribution of chemo-attractant density


The diffusion profiles of different emitters add up linearly. Therefore, the total distribution of chemo-attractant density from all present target cells can be written as

    F_2D(r⃗,t=t_n) = ∑_c=1^N_tarf_2D(|r⃗-r⃗_cn|).




 §.§ Modeling sensors for the chemo-attractant


In the case of temporal sensing, we assume that the immune cell can measure, in every time step n, the total density ρ^C_n = F_2D(r⃗=r⃗_C,t=t_n) of chemo-attractant at the center r⃗_C=r⃗_0,n of its cell body. It then computes the temporal difference

    Δρ^C_n = ρ^C_n - ρ^C_n-1.



In the case of spatial sensing, we assume that the immune cell has two sensors at the left and right border of its cell body, that is, at positions

    r⃗_L/R=r⃗_0,n+r_imm·(
    [ cos( ϕ_0,n-1±π/2 ); sin( ϕ_0,n-1±π/2 ) ]),



where the corresponding total chemo-attractant densities are ρ^L_n and ρ^R_n, respectively. It then computes the spatial difference 

    Δρ^LR_n = ρ^R_n - ρ^L_n.




 §.§ Mapping sensor signals to migration behavior


The two sensor signals available to the immune cell are the temporal density difference Δρ^C_n and the spatial density difference Δρ^LR_n.  The migration parameters which can be affected by these sensor signals are the speed of the immune cell v, its directional persistence ϵ, and its preference to turn left q_L.


For simplicity, we assume that the immune cell has two distinct migration modes, called the 'normal mode' N, and the 'approach mode' A. In the normal mode, the speed is v_N and the persistence is ϵ_N. In the approach mode, the speed is v_A and the persistence is ϵ_A. These four parameters can be tuned to optimize search performance.


At any time step n, the immune cell can only be in one of these two migration modes. The probability to be in the approach mode is computed as a function of the temporal gradient as follows

    q_A = prob(A) = logistic( c_A0+c_A1 Δρ^C_n ),

where

    logistic(x) = 1/1+e^-x

is the logistic function, and c_A0 as well as c_A1 are unknown coefficients that also have to be optimized. Note that for c_A1>0, the mode A is favoured whenever there is a positive temporal gradient, provided that the magnitude of the bias c_A0 is note too large.


In a similar way, the spatial gradient determines the probability q_R of the immune cell to turn right:

    q_R = prob(R) = logistic(c_R1 Δρ^LR_n ),

where c_R1 is an additional coefficient to be optimized. Note that for c_R1>0, right turns are favored whenever the chemo-attractant density at the right sensor is larger than that on the left sensor. 



 §.§ Choice of target cell density and linear system size


The density of target cells in the two-dimensional simulation plane is chosen to be

    ρ_tar := 1· 10^-5μ m^-2.

This density leads to a mean distance between nearest neighbors of

    r_nn = 1/2√(ρ_tar)≈ 158 μ m,

which is slightly larger than the diffusion length L_dif=100  μ m.


The linear system size is chosen as 

    L_sys = 1000 μ m,

which is considerably larger than L_dif and r_nn. The average number of target cells within the simulation area is

    N_tar = ρ_tar L^2_sys = 10.



Note that if an immune cell is migrating with its maximum speed of 6 μ m/min and with perfect directional persistence, it would take about 26 min ( = 26 simulation time steps) to cover the distance between two neighboring target cells. Within 100 min, an immune cell of perfect efficiency might encounter 3 to 4 target cells (ignoring the fact that r_nn is increasing slightly with each encounter and the simultaneous removal of the target). 



 §.§ Measuring search efficiency


We thus set the time period of a single simulation run to

    T_sim := 100   min.



After a specific simulation run k, the number of remaining target cells N_tar,k^rem is counted. We then quantify the efficiency of the immune cell by the number of eliminated target cells:

    Q_k = N_tar-N_tar,k^rem,

a quantity that can fluctuate considerable between each run. To overcome these fluctuations, the simulation is repeated

    N_runs := 10^4

times for each set of system parameters, using in each run a random initial configuration of the single immune cell and of the N_tar=10 target cells.


Finally, the search efficiency of the immune cell is defined as the average

    Q = 1/N_runs∑_k=1^N_runs Q_k.




 §.§ Optimization of response parameters



In general, the search efficiency in our model depends on up to n=7 unknown parameters π_i:

    Q = Q(π⃗) = Q(v_N, ϵ_N, v_A, ϵ_A, c_A0, c_A1, c_R1).



Finding the search strategy with the best search efficiency amounts to finding the parameter combination π⃗ that maximizes Q:

    π⃗_opt = arg max { Q(π⃗) }

We perform this quite high-dimensional numerical optimization using a grid-based variant of the 'Cyclic Coordinate Descent' method (CCD, see <cit.>). In each loop of this iterative method, the n parameters/coordinates π_i∈[1… n] are optimized one after the other in a cyclic way, greedily keeping the remaining n-1 coordinates at their presently best-performing values. In our variant of the method, an individual parameter π_k is optimized by evaluating Q=Q(π_k,{π_i≠ k}) for all discrete values of π_k on a regular grid within predefined minimum and maximum values, that is π_k∈[π_k,min, π_k,min+Δπ_k, … , π_k,max] The method stops when the same set of n optimal parameters is found in two subsequent iteration loops. 



 §.§ List of standard parameters


In Tab.(<ref>), we provide a list of all relevant system parameters, here called the Standard Parameters (SP). The first 12 parameters of the list are fixed for all simulations. During the optimization phase, a different random value of v_tar and ϵ_tar is drawn for each target cell, from uniform distributions in their respective ranges. During the evaluation phase, all N_tar target cells are set to the same values of v_tar and ϵ_tar, and these two parameters are then scanned through their ranges in subsequent simulation runs. The last 7 parameters of the list (v_N, ϵ_N, v_A, ϵ_A, c_A0, c_A1, c_R1) are free to be optimized within the given ranges.





§ RESULTS




 §.§ Concentration profile of chemo-attractant


In the fast diffusion limit, the global concentration distribution of chemo-attractant is a linear superposition of 'kernels', centered around the target cells. These kernels are the temporally stationary, rotationally symmetric solutions f_2D(r) of Eq.(<ref>). We have numerically computed the kernel for different diffusion lengths L_dif (See Fig.<ref>(b), in which the the green line corresponds to the case of standard parameters). As expected, the concentration profile decays almost exponentially for large radial distances r→∞.




 §.§ Blind search (BLS)


We start with an immune cell that completely lacks the ability to sense concentration gradients (c_A1=c_R1=0), and which is therefore performing a 'blind' search process (BLS). At the same time, we assume an extreme bias for the normal migration mode (c_A0=-5), which reduces the probability of the immune cell being spontaneously in the approach mode to an almost negligible value of q_A≈0.007 (The migration parameters of the approach mode are set to medium values v_A=3 and ϵ_A=0.5). The migration of such an immune cell can therefore be described as a homogeneous, correlated random walk with a fixed speed v_N and a fixed degree of directional persistence ϵ_N.


The N_tar=10 target cells, which are assigned random positions and migration directions before each simulation run, are assumed to form a widely mixed ensemble with respect to their migration parameters. For this purpose, at the beginning of each simulation run, we draw the speed and persistence parameters of each target cell independently from uniform distributions in the ranges v_tar∈[0,6] and ϵ_tar∈[0,1], respectively.


We first set the migration parameters of the immune cell to medium values v_N=3 and ϵ_N=0.5. In this case, the trajectory of the immune cell is not able to explore a significant part of the simulation area, even when the available time span is increased from the standard setting T_sim=100 min to T_sim=500 min (See Fig.<ref>(d). For a video, see <cit.>(UBS.mp4)). Repeating the simulation N_run=10^4 times, each spanning an evaluation period of T_sim=100, we find that the number Q_k of encounters between the immune cell and target cells is fluctuating from one run k to the next. The distribution p(Q_k) has an approximately exponential shape: In most simulation runs, the immune cell does not find any target, rarely one target, and almost never two targets. The average number of encounters with target cells, defined above as the search efficiency, is Q=0.110 in this case. If we let the immune cell migrate faster, using the parameters v_N=6 and ϵ_N=0.5, the search efficiency increases to Q=0.173. Additionally making the immune cell more directionally persistent, using the parameters v_N=6 and ϵ_N=1, results in a further increase of the search efficiency to Q=0.271. This demonstrates that even a blind, homogeneous search process can be optimized via the migration parameters v_N and ϵ_N. 


We therefore use CCD optimization to find the perfect migration parameters for the immune cell, again using the mixed ensemble of target cells throughout the optimization phase. It turns out that a blind, homogeneous search within a mixed ensemble of targets has the best efficiency Q when it is performed with maximum possible speed (in our case v_N=6) and with perfect directional persistence ϵ_N=1 (Fig.<ref>(a)).


The resulting optimal efficiency Q_BLS=0.27 can be seen as the overall performance of the immune cell, averaged over many possible types of target cells. In practice, it will also be of interest how the immune cell is performing against targets with specific, fixed  migration parameters. To investigate this 'versatility' of the immune cell, we have computed the search efficiency Q=Q(v_tar,ϵ_tar) of the optimized immune cell (that is, using v_N=6 and ϵ_N=1), as a function of the speed and persistence of the target cell (Fig.<ref>(b)). Here we find that the resulting search efficiency can vary between Q_min≈0.25 and Q_max≈0.32, depending on these two parameters. In particular, blind, homogeneous search works best when the targets are themselves fast and directionally persistent. Yet, if the targets exceed the immune cell with respect to the migration parameters, it is more appropriate to say that the targets are finding the immune cell than vice versa.


It is instructive to inspect the trajectory of the immune cell (Small gray dots in Fig.<ref>(c). For a video, see <cit.>), in relation to the targets, over an extended time period. For this purpose, we set the speed of the targets to zero, so that they remain stationary throughout the entire simulation. Since the persistence of the optimized immune cell is ϵ_N=1 in the normal mode, the trajectory is straight for most of the time (Note the effect of periodic boundary conditions). However, with a tiny probability of q_A≈0.007, the immune cell also adopts the 'approach mode', where the migration parameters are v_A=3 and ϵ_A=0.5, and these rare events lead to an abrupt change of direction. It is remarkable that there occur several 'near misses' between the immune cell and one of the targets. Yet, without any sensing abilities, the immune cell most of the time cannot seize these opportunities.



 §.§ Random mode switching (RMS)


We continue to consider blind search, characterized by the absence of sensitivity for concentration gradients  (c_A1=c_R1=0). But this time we allow the immune cell to switch between its two migration modes randomly and spontaneously, a situation that creates a heterogeneous correlated random walk. For this purpose, we now declare not only the parameters v_N and ϵ_N, but also c_A0, v_A and ϵ_A as free, optimizable parameters.


Although the system is now considerably more flexible than in the case of homogeneous blind search, CCD optimization shows that this flexibility brings no significant improvement of the the search efficiency (Fig.<ref>(d)), as Q_RMS=0.28≈0.27= Q_BLS. Indeed, the optimal efficiency is found for a bias c_A0=5, which keeps the immune cell in the approach mode virtually all the time, thus leaving the values v_N=4 and ϵ_N=0 irrelevant. Within the approach mode, the optimized immune cell is as fast (v_A=6) and persistent (ϵ_A=1) as possible, just like in the above  homogeneous BLS case. This demonstrates that in blind search, purely spontaneous mode switching performs worse than a homogeneous random walk at maximum speed and perfect directional persistence. Since the optimal RMS strategy is - except for a name change of the dominating migration mode - identical to the BLS strategy, we also find the same results for Q=Q(v_tar,ϵ_tar) (Fig.<ref>(e)). The sample trajectory of the immune cell also resembles that of the BLS strategy (Fig.<ref>(f). For a video, see <cit.>).



 §.§ Temporal gradient sensing (TGS)


Next, we investigate how the killing efficiency can be enhanced when the immune cell is able to measure temporal gradients of the chemo-attractant and to switch between the normal mode N and the approach mode A accordingly. In order to make this adaptive mechanism work, there are six parameters to be optimized: The speeds (v_N,v_A) and persistence values (ϵ_N,ϵ_A) in the two migration modes, as well as the bias of the approach mode (c_A0) and the sensitivity for temporal chemo-attractant gradients (c_A1). Without restriction of generality, the latter quantity is assumed to be non-negative, c_A1≥ 0, because a positive temporal gradient of the chemo-attractant Δρ^C_n means that the immune cell is approaching a target cell, and this should increase the probability of the approach mode A, whatever this means for the speed and persistence of the immune cell.


CCD optimization shows (Fig.<ref>(a)) that the optimum bias for the approach mode is c_A0=2, which corresponds to a probability q_A≈0.88 of the immune cell being in the approach mode if detecting no or only a very weak temporal gradient. When however a significant gradient is present, the large sensitivity parameter c_A1=500 causes an almost deterministic mode switching behavior: In positive gradients, the optimal TGS cell is adopting the approach mode, which is maximally fast (v_A=6) and persistent (ϵ_A=1). In negative gradients, it is adopting the normal mode, which is also fast (v_N=6), but directionally non-persistent (ϵ_N=0). The resulting search efficiency of the optimized TGS strategy against target cells with mixed migration properties is Q_TGS=1.07, which surpasses the blind strategies by a factor of Q_TGS/Q_BLS≈ 4. 


Confronted with target cells of fixed migration properties (Fig.<ref>(b)), the performance of the optimized TGS strategy is degrading relatively quickly when the targets are fast and directionally persistent. 


The sample trajectory of the immune cell (Fig.<ref>(c). For a video, see <cit.>) demonstrates the alternating phases of zero persistence (ϵ_N=0, 'zigzag'-like motion) and perfect persistence (ϵ_A=1, straight motion). In contrast to the blind search strategies, the immune cell is now able to perfectly home in on a target, once it came close to it.  



 §.§ Spatial gradient sensing (SGS)


We now consider an immune cell that is virtually always in the approach mode (enforced by c_A0=500), but has the ability to turn left (clock-wise) or right in response to the spatial gradient of the chemo-attractant. The relevant response coefficient for this mechanism is the sensitivity c_R1. Yet, how well the immune cell can follow a spatial gradient also depends on the migration parameters v_A and ϵ_A, because they determine how quickly the cell can adjust its direction as it follows a spatial gradient. 


CCD optimization shows (Fig.<ref>(d)) that the optimized SGS immune cell turns into the direction of larger chemo-attractant concentration with maximum sensitivity (C_R1=500). It migrates with maximal speed (v_A=6), but with a specific degree of persistence that is smaller than one (ϵ_A=0.8).  The resulting search efficiency of the optimized SGS strategy against target cells with mixed migration properties is Q_SGS=2.58, which surpasses the TGS strategy by a factor of  Q_SGS/Q_TGS≈ 2.4, and blind strategies by a factor of Q_SGS/Q_BLS≈ 9.6.


Confronted with target cells of fixed migration properties (Fig.<ref>(e)), the optimized SGS cell has a relatively constant performance for targets with small to medium speeds and persistences. In the extreme case of targets with v_tar≈ 6 and ϵ_tar≈ 1, the performance declines, but even then it is still about as good as the optimal TGS performance. 


The sample trajectory (Fig.<ref>(f). For a video, see <cit.>) shows that the optimized SGS immune cell is wasting almost no time between subsequent target attacks. It moves from one target to the next in an efficient way, resembling the optimal solutions of a traveling salesman problem.



 §.§ Combined spatial and temporal gradient sensing (CGS)


Finally, we consider an immune cell that can, both, switch between two migration modes in response to the temporal chemo-attractant gradient, and at the same time turn left and right in response to the spatial gradient. Since these two mechanisms have different requirements with respect to the migration parameters (For example, temporal sensing requires ϵ_N=0, but spatial sensing works best with ϵ_N=0.8), it is not clear whether a combination of the two abilities is advantageous or reduces the killing efficiency.


CCD optimization of combined gradient sensing involves the complete set of eight free parameters (Fig.<ref>(g)). The resulting bias c_A0=5 means that the optimized CGS cell is adopting the approach mode practically all the time. In this mode, it just performs spatial gradient sensing, since all the parameters that are relevant to SGS are actually unchanged (c_R1=500, v_A=6, and ϵ_A=0.8). However, the optimized CGS cell is also highly sensitive to temporal gradients (c_A1=500). Therefore, in the presence of a sufficiently negative temporal gradient, it will switch to the normal mode, which is fast (v_N=6) but only medium persistent (ϵ_N=0.5). This means that Combined gradient sensing is basically like spatial gradient sensing, but with the additional feature of a less persistent migration in strongly negative temporal gradients. The resulting search efficiency of the optimized CGS strategy against target cells with mixed migration properties is Q_CGS=2.61, which is only slightly better than the SGS strategy. The versatility of combined gradient sensing resembles that of purely spatial gradient sensing (Fig.<ref>(h)). Also the sample trajectory (Fig.<ref>(i). For a video, see <cit.>) has basically the same characteristics as in the SGS strategy.



§ DISCUSSION AND SUMMARY


Chemotaxis is an important phenomenon in prokaryotic cells, such as bacteria <cit.>, but it plays an equally fundamental role in eukaryotes <cit.>. Providing one of the simplest examples of goal-directed behavior, chemotaxis is also a fascinating topic that has attracted the interest of physicists since decades.

Among the physics-oriented publications, many focus on the formation of spatio-temporal patterns in colonies of self-propelled agents with mutual predator-prey relations <cit.>. While most of these studies describe the complete density distribution of the agents by partial differential equations, as in the Keller-Segel model, a few works also describe the motion of the agents individually, for example in the framework of active Brownian particles <cit.>. A more recent theoretical study <cit.> investigates the chemotactic pursuit of a single prey agent by a predator. Although this work addresses a research question similar to ours, it is based on different model assumptions. In particular, it assumes not only that the predator is chemically attracted by the prey, but also that the prey is repelled from the predator. Furthermore, the guiding chemicals in <cit.> are assumed to have an infinite life time, which prevents the formation of a stable chemical 'cloud' around each agent and leads to long-range interactions.

Most existing physics-oriented models of single-agent chemotaxis describe the motion of the predator by a Langevin equation of an over-damped particle that feels a delta-correlated random force, plus a deterministic force proportional to the local concentration gradient of the guiding substance. According to this model, in a situation with zero gradient, the predator will perform a Gaussian random walk without any directional correlations. In the presence of a  non-zero gradient, a force will act on the predator that pulls it into the direction of the prey.

In reality, however, cells typically move with a high degree of directional persistence, even without detecting any chemical concentration differences. Also, changing the migration direction of a cell is a complex process that involves, among many other stochastic events, a partial remodeling of the cytoskeleton. Even though a concentration gradient is pointing strongly into a particular direction, the predator may not be able to turn into this new direction immediately and with sufficient precision. Instead of describing chemotaxis as a deterministic force, an indirect, stochastic approach may therefore be more appropriate: The migration of the predator is a stochastic process, and the parameters of this process (such as the probability of turning left or right) are modified as a function of the chemotactic gradient.

Apart from this opportunity to improve the modeling of chemotaxis, the present work was motivated by preliminary experiments <cit.> on chemotactic 'pursuit' in a Petri dish. These experiments studied the interaction of natural killer (NK) cells <cit.>, extracted from the blood of human donors, and human leukemic K562 tumor cells. By following the migration path of all cells over several hours, it was observed that a fraction of the highly mobile NK cells approached individual tumor cells (which themselves showed only weak mobility) in a directionally persistent walk, attacked them, and often induced their death subsequently. These in-vitro experiments demonstrated that single immune cells are able to find certain pathogenic targets on their own account. However, it remained unclear if the observed attacks were merely chance encounters, or actually guided by chemotactic mechanisms.

In order to answer this question, it would certainly be useful to know which efficiencies can be expected from a 'blind' search, and from different 'guided' search strategies based on chemotaxis. In addition, finding distinct, highly efficient and robust search strategies by numerical optimization of a simulated immune cell would also reveal characteristic patterns of search behavior, that might then be used as characteristic 'fingerprints' of goal-directed search in future automatic detection algorithms. 

For these reasons, we have compared in this work five distinct strategies of search, namely blind search with fixed speed and directional persistence of the immune cells (BLS), blind search with random switching between two distinct migration modes (RMS), guided search based on temporal gradients of the chemo-attractant (TGS), guided search based on spatial gradients of the chemo-attractant (SGS), and a combination of temporal and spatial sensing (CGS). Throughout our study, we have kept the system geometry (two dimensions, as on a Petri dish) and all parameters (density of the target cells, properties of the chemo-attractant, sizes and migration properties of the cells, sensing abilities of the immune cell) close to experimentally realistic values.

In the case of blind search (BLS), not surprisingly, the search efficiency of the immune is almost an order of magnitude lower than with the best guided search mechanisms. Nevertheless, since many pathogens will not emit any chemical substance that the immune cell can detect and use as a guide to its target, blind search may often be the only option. It is therefore fortunate that blind search can be easily optimized by making the immune cell as fast and directionally persistent as possible. This can be understood most easily assuming immobile target cells that are located at random positions within the plane. As the search time t is going on, the blindly migrating immune cell is exploring more and more regions of the Petri dish, and we can mentally mark all spatial pixels that have been visited at least once by the immune cell. The total area of all marked pixels, A(t), here called the 'visited area', is growing monotonously with time, and all target cells that happen to be located within the visited area can be considered as found by the immune cell. Their expected number is N_found(t)=ρ_tar A(t), where ρ_tar is the areal density of target cells. If the immune cell is migrating with low directional persistence, it will re-visit many pixels more often than once, which is counter-productive with target cells that never move. In this case, the visited area will grow sub-linearly with time. By contrast, A(t)∝ t for an immune cell that is migrating with perfect directional persistence and constant speed, that is, uniformly along a straight line. It is therefore clear that high directional persistence is an important way to improve the blind search efficiency Q of immune cells. At the same time, speed is another key factor for efficient search: For an immune cell in uniform motion, the expected number of found target cells at the end of the search period, N_found(t=T_sim)∝ v_imm, will be directly proportional to its speed.

However, it is well-known that actual cells - and not only immune cells - are showing gradual or abrupt changes of their speed and persistence <cit.>, so that their migration has to be described by a temporally heterogeneous stochastic process. The result of such parameter fluctuations are often 'anomalous' properties of the cell's random walk, such as a mean squared displacement that increases with lag-time approximately as a powerlaw. It is not clear whether  temporally heterogeneous cell migration is just a side effect of other causes (such as differences in the local micro-environment of the migrating cell or internal changes connected with the cell cycle), or if it actually serves a purpose. Theoretically, the heterogeneity may help to increase the blind search efficiency of an immune cell, particularly when the targets are mobile. We have therefore investigated how the search efficiency is affected when the immune cell performs random switches between two different migration modes (the RMS strategy). Yet, as suggested by the theoretical argument above, the RMS strategy did not perform significantly better than blind search with fixed migration parameters.

Next, we have investigated guided search strategies that are based on the sensing of chemotactic gradients. In the case of temporal gradient sensing (TGS), we found that the optimized immune cell is switching between two distinct migration modes: In positive gradients, it is adopting the approach mode, which is maximally fast (v_A=6) and persistent (ϵ_A=1). By this way, the cell is climbing up the gradient consistently, which usually corresponds to approaching one of the targets. In negative gradients, it is adopting the normal mode, which is also fast (v_N=6), but directionally non-persistent (ϵ_N=0). In this mode, the cell is exploring new migration directions, until it finds one with a positive gradient. Note that the optimal TGS strategy found here by numerical parameter optimization strongly resembles the chemotaxis behavior of Escherichia Coli <cit.>, with its gradient-dependent switching between swimming and tumbling modes of migration. Compared to blind search, TGS is more effective. On the other hand, the gained factor of four in search efficiency is not really large.

In the case of spatial gradient sensing (SGS), we found that the optimized immune cell turns into the direction of larger chemo-attractant concentration with maximum sensitivity (C_R1=500). It migrates with maximal speed (v_A=6), but with a specific degree of persistence that is smaller than one (ϵ_A=0.8). Presumably, this specific degree of persistence represents an optimal compromise between the need to maximize the visited area, and the need to perform clockwise and counter-clockwise turns with the right curvature. Compared to blind search, SGS is almost an order of magnitude more efficient. A combination of temporal and spatial sensing (CGS) turned out to bring no significant advantages compared to pure spatial sensing. 

The blind and guided search strategies differ characteristically in how the search efficiency Q=Q(v_tar,ϵ_tar) depends on the migration parameters of the targets: While blind search (BLS, RMS) works better with fast and persistent targets, the opposite is true for guided search (TGS, SGS, CGS). In guided search, due to the optimization against a mixed set of targets, the search efficiency Q=Q(v_tar,ϵ_tar) remains approximately constant for most combinations of v_tar and ϵ_tar. Only for targets that are simultaneously extremely fast and persistent does Q decline significantly. Assuming an experimental possibility to vary the migration properties of the targets, without affecting the immune cell or the properties of the chemo-attractant, this predicted difference in Q=Q(v_tar,ϵ_tar) offers a first possibility to distinguish between blind and guided search strategies.

Finally, our work suggests how to detect different search strategies of an immune cell by looking for characteristic patterns in the cell's trajectory: In the case of temporal sensing, the immune cell will show alternating phases of low and high directional persistence, and the probability of the high persistence mode will increase whenever the immune cell approaches one of the targets. In the case of spatial sensing, the left- and right-turns of the immune cell will occur in such a way that they tend to align the cell into the direction of the closest target.

In future work, our investigation could be improved and extended in several obvious ways. For example, we have so far assumed that the immune cell is able to detect arbitrarily small concentrations (or differences between two concentrations) of the chemo-attractant. A lower detection limit may very well change the optimal search parameters and, accordingly, the associated search strategies. Furthermore, it would be straight-forward to extent our simulations to three spatial dimensions. Since the local concentration is, at least within the fast diffusion limit, just a sum over fixed kernels, a change from 2D to 3D kernels would not increase the computation time. By contrast, going beyond the fast diffusion limit is computationally more demanding, as it requires to solve the partial differential equation of the spreading and decaying chemo-attractant along with the motion of the cells. However, we have already demonstrated the feasibility of this approach in 2D (data not shown).  






§ ACKNOWLEDGEMENTS


We thank Christoph Mark and Franziska Hoersch for many helpful discussions, and for sharing preliminary experimental results on the interaction between natural killer cells and tumor cells in a Petri dish.



§ FUNDING

This work was funded by the Grant ME1260/11-1 of the German Research Foundation DFG.



§ REFERENCES


unsrt















