Detecting long-range interactions between migrating cells
Claus Metzner, Franziska Hörsch, Christoph Mark, Tina Czerwinski, Alexander Winterl, Caroline Voskens, Ben Fabry
DDetecting long-range interactions between migrating cells
C. Metzner * , F. H¨orsch * , C. Mark * , C. Voskens + , and B. Fabry ** Biophysics, Friedrich-Alexander University Erlangen-N¨urnberg, Erlangen + Department of Dermatology, University Hospital Erlangen, Friedrich-AlexanderUniversity Erlangen-N¨urnberg, ErlangenMarch 20, 2020
Abstract
Chemotaxis enables cells to systematically approach distant targets that emit a diffusibleguiding substance. However, the visual observation of an encounter between a cell and a targetdoes not necessarily indicate the presence of a chemotactic approach mechanism, as even ablindly migrating cell can come across a target by chance. To distinguish between chemotacticapproach and blind migration, we present an objective method that is based on the analysisof time-lapse recorded cell migration trajectories: For each movement of a cell relative to theposition of potential targets, we compute a p -value that quantifies how likely the direction ofthis movement step is under the null-hypothesis of blind migration. The resulting distribution of p -values, pooled over all recorded cell trajectories, is then compared to an ensemble of referencedistributions in which the positions of targets are randomized. We validate our method withsimulated data and demonstrate that the method reliably detects the presence or absence oflong-range cell-cell interactions. We then apply the method to data from experiments withhighly migratory human-derived natural killer (NK) cells and nearly stationary K562 tumorcells as targets that are dispersed in a 3-dimensional collagen matrix. Although NK cells findand attack K562 cells with high efficiency, we show that the attacks are consistent with atarget-blind random walk. NK cells achieve their high attack efficiency by combining highmigration speeds with high degrees of directional persistence. We provide a freely availablePython implementation of the p-value method that may serve as a general tool for extractingthe interaction rules in collective systems of self-driven agents. a r X i v : . [ q - b i o . Q M ] M a r ntroduction Pursuit and evasion are ubiquitous in nature [1]. An obvious example are predator-prey relations,where an agent attempts to catch a target that is struggling to escape. In such extreme cases, thecontinuous mutual reaction of the two opponents proves without doubt that they act in a goal-directed way and are able to sense each other from a distance. By contrast, when mobile agents areforaging for non-evading targets, it is not always clear whether the agents are performing a blindrandom walk and find their targets merely by chance, or if they recognize them already from largerdistances and approach them systematically. The situation is particularly ambiguous in the case ofmicro-organisms, which can only migrate with a relatively large degree of directional randomness [2],and for which chemotaxis is often the only available mechanism to locate distant targets [3].We have encountered such ambiguous behavior in experiments with highly mobile natural killer (NK)cells and almost immobile tumor cells, randomly distributed inside a 3-dimensional collagen gel (fordetails see below). By continuously monitoring the migration paths of the cells, we regularly observeNK cells that migrate to a nearby tumor cell, establish steric contact and attack the tumor cell,causing its subsequent death. This killing behavior after establishing cell-cell contact is consistentwith the expected function of the immune system to eliminate pathogens. However, the interpretationof the preceding phase, in which the NK cells approach their targets, is ambiguous: On the one hand,it is known that NK cells can be chemotactically recruited by other cells of the immune system thathave located a pathogen [4, 5]. Given this fundamental chemotactic ability, the NK cells mightvery well be able to follow chemical traces of tumor cells directly. On the other hand, thoroughvisual inspection of time-lapse video recordings yields numerous examples where a migrating NKcells misses or even seems to turn away from a nearby tumor cell.The question of how immune cells are locating pathogens is of general importance, e.g. for optimizingcell-based immunotherapies for cancer [6, 7, 8]. A method for quantitatively analyzing the strengthand range of chemotactic interaction could ultimately help opening ways to modify and improve theforaging efficiency of the immune cells by external interventions.In this study, we present a method for detecting long-range interactions between two types of agents,based only on their time-lapse recorded trajectories. The method could, in principle, be applied toany kind of self-propelled agents, including multi-cellular organisms or animals. However, we focushere on unicellular organisms, and specifically on immune cells and their targets. Our approach is asfollows: We view the detection of interactions as a statistical test of the hypothesis that immune cellmigration is affected by distant targets , against the null-hypothesis that the immune cells performa blind random walk . Every discrete step of an immune cell’s migration trajectory is characterizedby a p-value, describing the probability that the step is part of a blind random walk. Computingthe distribution of p-values over all recorded immune cells and time steps then reveals the presenceof long-range interactions between immune cells and targets by a statistically significantly largerfraction of small p-values compared to a reference distribution of p-values in which the positions ofthe targets are randomized.We validate our method with surrogate data, using a simulation framework for chemotactical behaviorthat has been published previously [9]. In a first simulated scenario, called ’blind search’ (BLS), themigration of immune cells is not influenced by the presence or the positions of targets. The immunecells migrate blindly with individually different but temporally constant migration parameters. Inthe second case, called ’random mode switching’ (RMS), immune cells are also migrating blindly, buttheir migration behavior occasionally switches e.g. from slow to fast or from random to persistent,for reasons unrelated to the targets - a scenario that has been previously shown to be ubiquitous[10]. In the third test case, called ’temporal gradient search’ (TGS), immune cells are able to detectdifferences of the target-related chemo-attractant concentration over time, and they modulate theirdegree of directional persistence accordingly - a well-known strategy of chemotaxis that is found2n the run-and-tumble behavior of E.coli [11]. In the fourth test case of ’spatial gradient search’(SGS), immune cells can directly measure and turn in in the direction of a spatial gradient in theconcentration of a chemo-attractant. Applied to these four simulated scenarios, the p-value methodcorrectly finds the absence of long-range interactions between immune and target cells in the cases of’blind search’ and ’random mode switching’, and the presence of interactions in the cases of ’temporalgradient search’ and ’spatial gradient search’.Moreover, we apply the p-value method to data from experiments with human natural killer (NK)immune cells and K562 tumor cells inside a three-dimensional collagen gel. NK cells are able to findand kill a substantial fraction of the tumor cells over a time course of several hours. Despite this highkilling efficiency, the resulting p-value distributions are consistent with the null-hypothesis that NKcells perform a target-blind random walk. We explain this apparent contradiction with the specificmigration properties of NK cells: Due to their high speed of typically 6 µ m/min [12], they can covera larger search area than most other cell types in a given time period. Moreover, their relatively highdirectional persistence makes repeated visits of the same locations less likely, which increases thesearch efficiency [13]. We speculate that a fast and directionally persistent migration behavior hasevolved in immune cells as a generally effective, target-independent search strategy against a varietyof possible target types, including those that do not emit any detectable chemo-attractants. Materials and Methods
Step 1: Data generation
Experimental setup
We generally assume an experimental assay where immune and cancer cells are mixed together in a3-dimensional matrix that is suitable for effective cell migration and enables proper imaging with amicroscope. If the matrix layer has a vertical thickness of only a few cell diameters, the system canbe considered quasi two-dimensional, and the subsequent analysis can be restricted to the horizontal(x,y) cell positions. In the case of thicker matrices, the z-positions of the cells need to be measured aswell to select and analyze only those pairs of immune and target cells that are approximately in thesame horizontal plane. In principle, our method can be extended to include also the z-coordinates ofcell migration trajectories, but because of the non-isotropic resolution of most microscopes, we limithere our analysis to the horizontal plane. We assume that the microscope’s field of view is time-lapserecorded with a constant time interval ∆ t between successive recordings. This ∆ t has to be shortenough ( ≈ ≈
500 frames), and with large numbers of cells in the field of view ( ≈ · NK immune cells and 3 · K562 tumor cells are mixed withice-cold 500 µ l acid-dissolved collagen solution (1.2mg/ml) and pipetted into a tissue-culture-treated35 mm dish (for a detailed protocol, see [10]). The polymerization of the collagen solution is initiatedby placing the dish for 30 min in a cell culture incubator at 37 ◦ C, 5% CO . Due to surface tension,the thickness of the polymerized collagen gel decreases towards the center of the dish with a heightof ≈ µ m (Fig. 3(j,k)). Time-lapse imaging can thus be realized in bright-field mode withoutscanning in z-direction, while the cells still showed the same characteristic migration behavior as ina thick collagen gel. We recorded 9 independent data sets, each including between 333 and 1547images, with a time interval of 45 seconds between two subsequent frames. The images had 1344 x1024 pixels with a linear size of 0.645 µ m. 3 ell tracking To extract the information required for our interaction-detection algorithm, the recorded cells need tobe detected and individually tracked, yielding the 3D center-of-mass coordinates (cid:126)R ( i ) t = ( x ( i ) t , y ( i ) t , z ( i ) t )of every cell i in each video frame t . Each cell needs to be labeled with a unique ID number i thatremains consistent over subsequent frames. Furthermore, each cell must be classified as either c = 0(immune cell) or c = 1 (target cell). All information regarding a particular cell at a particulartime needs to be stored in an ’observation’ vector of the form ( t, x, y, z, i, c ). The total numberof observations may change between frames, as cells may leave or enter the microscope’s field ofview, because of cell division and death, or also due to tracking problems. The observations fromall video frames (in any order) need to be combined into a matrix, with each row corresponding toan observation vector. This matrix, stored as a Numpy-array, forms the input to our interaction-detection algorithm.In this study, cells are automatically segmented using local differences of image entropy (Fig. 3(l)).The classification in immune and tumor cells is based on differences in speed, size, and brightness.Finally, the temporal trajectory of each cell is determined using the overlap of the cell area betweensuccessive time frames. Generation of simulated data
To validate our p-value method of interaction detection, we use a previously published softwareframework for the simulation of chemotactic hunting behavior [9], which provides the following fourscenarios: (1) In ’Blind Search’ (BLS), the immune cells do not interact at all with the targets butmigrate blindly, according to a correlated random walk with fixed parameters for the mean step width(speed) and for the degree of directional persistence. (2) In ’Random Mode Switching’ (RMS), theimmune cells are still blind with respect to the targets, but occasionally switch between a highly per-sistent random walk and a non-persistent (diffusive) random walk mode. (3) In ’Temporal GradientSensing’ (TGS), the immune cells actually approach the targets by following the temporal gradientof chemo-attractant. In particular, the model assumes that the immune cells perform a highly per-sistent correlated random walk as long as the concentration of chemo-attractant is increasing withtime. When the concentration is decreasing, the immune cells switch to a diffusive (uncorrelated)random walk in order to find a more goal-directed migration direction. (4) In ’Spatial GradientSensing’ (SGS), the immune cells are able to measure the spatial gradient of chemo-attractant andto actively turn into the direction of a nearby target.To generate the surrogate data for the present paper, we set all parameters of the chemotacticsimulating framework to the same values as in [9]. However, while all results were averaged over10000 runs in [9], we now produce for each scenario only a single simulated data set with a longerduration (500 time steps of ∆ t = 1min), with a larger field of view (5000 µ m x 5000 µ m) and witha larger number of cells (100 immune cells and, initially, 50 target cells). Step 2: Data filtering ’Triplet’-based analysis
The elementary unit for our data analysis is a ’triplet’, consisting of three consecutive cell positions (cid:110) (cid:126)R ( i ) t − , (cid:126)R ( i ) t , (cid:126)R ( i ) t +1 (cid:111) of an immune cell i (solid blue dots in Fig. 1(a)). Our algorithm automatically ex-tracts all such triplets from the matrix of observations. Recording or tracking gaps are automaticallyexcluded from the list of triplets. 4 xcluding cells with too few triplets From the list of all triplets, we extract the subset belonging to a particular immune cell i . If this cellhas a number N ( tr ) i of triplets smaller than N ( tr ) min (typically set to 5), the cell is excluded from thesubsequent analysis, because it is not possible to reliably estimate the average migration behavior ofa cell based on such a small number of positions (for details see below). Excluding triplets with too distant targets
Often, there will be a-priori knowledge about the maximum expected distance r max for interactionsbetween immune and target cells (Otherwise, r max can be set to a value larger than the size of themicroscope’s field of view). For the purpose of interaction detection (for details see below), we restrictthe analysis to triplets which have at least one target present in a sphere of radius r max around thetriplet’s central immune cell position (cid:126)R ( i ) t = ( x ( i ) t , y ( i ) t , z ( i ) t ). Generally, a smaller r max reduces thecomputation time of the algorithm, but a larger r max increases the number of possible targets andthus reduces statistical fluctuations in the evaluation. Step 3: Data analysis
Cell migration model
We focus on the in-plane, horizontal motion of the cells. For this purpose, we only use two-dimensionalcoordinates, in the following denoted by lower-case position vectors (cid:126)r ( i ) t = ( x ( i ) t , y ( i ) t ). The sequenceof horizontal positions (cid:126)r of each individual cell i over successive time indices t = 0 , , , . . . is ap-proximated by a directionally persistent, discrete time random walk. It is characterized by a certaindistribution p i ( w ) of step widths w , and a distribution p i ( θ ) of turning angles θ . Here, the step width(Euclidean distance) for a cell’s movement between time t and t +1 is defined as w = | (cid:126)r ( i ) t +1 − (cid:126)r ( i ) t | , andthe turning angle is defined as the angle between the two shift vectors (cid:104) (cid:126)r ( i ) t +1 − (cid:126)r ( i ) t (cid:105) and (cid:104) (cid:126)r ( i ) t − (cid:126)r ( i ) t − (cid:105) .The step width distribution is modeled as a Rayleigh distribution with speed parameter σ i : p i ( w ) = wσ i exp (cid:18) − w σ i (cid:19) . (1)The turning angle distribution is modeled as a von-Mises distribution with persistence parameter κ i : p i ( θ ) = 12 πI ( κ i ) exp ( κ i · cos( θ )) . (2)We estimate the individual speed parameter σ i and persistence parameter κ i for each cell i , basedon its complete recorded time series, as described in [14]. These two parameters characterize theaverage in-plane migration behavior of the cell. ’Ordinary’ and ’extraordinary’ steps In our method of interaction detection, we focus on the turning angles θ of the immune cells, aquantity that is statistically fluctuating from one step to the next, approximately described by thevon-Mises distribution with persistence parameter κ . For positive κ , corresponding to directionallypersistent migration, the von-Mises distribution is peaked around a zero turning angle, so that mostof the turning angles will have small magnitudes (’ordinary moves’) and only relatively few will havelarge magnitudes (’extraordinary moves’).For the sole purpose of visualization, we can set an arbitrary threshold angle θ thr and define ordinarymoves as those with | θ | ≤ θ thr . For example, the threshold θ thr could be chosen such that ordinarymoves occur in target-blind migration with a probability of P ord = 0 .
95. Graphically, the interval ofordinary turning angles can then be depicted as a ’persistence cone’ (Fig. 1(a,b)). Turning angles5hat lead to the outside of the persistence cone would then be regarded as extraordinary (Case (3) inFig. 1(b)). Note, however, that our method of interaction detection is directly based on the von-Misesdistribution, and neither θ thr nor P ord play any role in the calculation of the p-values. The definitionof the persistence cone is only used to illustrate the fundamental idea of the method. Evidence for target pursuit
Even if an immune cell is moving with high persistence into the direction of a nearby target (reddots in Fig. 1(b)), this provides no evidence for target pursuit if the individual migration steps areclassified as ordinary (cases (1) and (2) in Fig. 1(b)). Only steps that are classified as extraordinaryand at the same time are highly target-directed provide some evidence for target pursuit (case (3) inFig. 1(b)), in particular if they occur more often than would be expected for a target-blind immunecell.
Definition of the ’approach cone’
To quantify the target-directedness of a step, we define an ’approach cone’ (left gray area in Fig. 1(a))as follows: We consider three successive positions of an immune cell i (solid blue circles), given bythe triplet (cid:110) (cid:126)r ( i ) t − , (cid:126)r ( i ) t , (cid:126)r ( i ) t +1 (cid:111) , as well as the position (cid:126)r ( j ) t of a nearby target cell j (solid red circle).The optimal turning angle θ (cid:63) (dashed red arc) of a goal-directed immune cell shifts the migrationtrajectory directly towards the target cell with an optimal shift vector (cid:126)s opt = (cid:126)r ( j ) t − (cid:126)r ( i ) t (red vector). Inpractice, immune cell i has moved along the vector (cid:126)s = (cid:126)r ( i ) t +1 − (cid:126)r ( i ) t (solid blue vector), which enclosesan angle ∆ θ (red arc) with the optimal shift vector (cid:126)s opt . There exists another, hypothetical shiftvector (cid:126)s that encloses the same angle ∆ θ with the optimal shift vector (cid:126)s opt (dashed blue vector).The interval of turning angles between θ = θ (cid:63) − ∆ θ and θ = θ (cid:63) +∆ θ ia the approach cone, that is,the set of directions which are at least as target-oriented as the actual shift of the immune cell. Definition and interpretation of the p-value
By integrating the von-Mises distribution p i ( θ ) over all turning angles θ ∈ [ θ (cid:63) − ∆ θ, θ (cid:63) +∆ θ ] withinthe approach cone, we compute a p-value, subsequently denoted by the symbol ˆ p (green area underthe p ( θ ) curve in Fig. 1(a)). ˆ p = (cid:90) θ (cid:63) +∆ θθ (cid:63) − ∆ θ p i ( θ ) dθ. (3)ˆ p can be interpreted as the probability that the observed move of the immune cell, or an even moretarget-directed move, could occur in a target-blind migration . Very small p-values indicate thatimmune cells are attracted towards target cells, while very large p-values indicate that immunecells are repelled from target cells. Due to its definition, the p-value can be very small only ifthree conditions are simultaneously fulfilled: (A) The persistence cone is narrow (high directionalpersistence of the immune cell, corresponding to a narrow von-Mises distribution p i ( θ )). (B) Theapproach cone is narrow (the immune cell turns almost exactly towards the target cell). (C) Thetwo cones are non-overlapping and distant from each other (the immune cell is ’going out of its way’to approach the target). Distribution of p-values
Finding just a few steps with very low p-value does in general not provide convincing evidence fora target-directed immune cell migration. Moreover, a subset of immune cells might be attractedto the targets while others are repelled from them. Alternatively, the same immune cell could beattracted and repelled by targets at different times. All these cases are comprehensively described bythe global probability distribution q obs (ˆ p ) of all observed p-values. We approximate this continuousprobability distribution by a discrete histogram (Fig. 2 and Fig. 3).6 eference distributions of p-values Since extraordinary steps occur also in target-blind migration with probability 1 − P ord , and somesome of these steps may accidentally lead into the direction of nearby targets, we need to compare q obs (ˆ p ) with a reference distribution q ref (ˆ p ) of a system that resembles the observed one in all respects,except that there are no interactions between immune cells and target cells. To obtain this referencedistribution, we use the following bootstrapping method [15]: We generate a reference data set bykeeping the positions of the immune cells unchanged but shifting all target cells that are locatedwithin the maximum interaction radius r max to new, independent random positions within thatradius (Fig. 1c). From this reference data set, we compute the histogram of p-values, yielding a(first) reference distribution q ref (ˆ p ) with the same sample size as q obs (ˆ p ). Since we are interestedin the fluctuations of the q ref (ˆ p )-values in each histogram bin, we repeat the same procedure for alarge number N s ≈
100 of statistically independent reference data sets (Thin gray lines in Fig. 2and Fig. 3). From these N s histograms, we compute the mean µ k and standard deviation σ k ofthe q ref (ˆ p )-values in each histogram bin k . Based on this statistics, we define confidence intervals(dashed blue and red lines in Fig. 2 and Fig. 3) for each bin k as [ µ k − . σ k , µ k + 1 . σ k ].Assuming a normal distribution, the probability of a value above the upper (or below the lower) limitof the confidence interval is then 0.05 in the target-randomized reference system. If the measuredp-value distribution q obs (ˆ p ) lies outside the confidence interval of the reference systems at least insome histogram bins, this may be interpreted as a statistically significant effect, indicating that thetargets somehow affect the migration of the immune cells.In our case, the width ∆ q k of the confidence intervals has been arbitrarily set to 2 × . σ k , becausethe resulting significance level of 0.05 is a common choice in science. The user is however free tochoose other values for ∆ q k , such as ∆ q k = 2 × . σ k for a significance level of 0.01. Results
Validation of p-value method with simulated data
We first validate the p-value method with data from chemotaxis simulations, using a recording periodof ∆ t = 1 min, a maximum detection radius of r max = 500 µ m, and a number of N s = 100 referencedistributions. In the BLS scenario, as expected, the resulting p-value distribution is almost identicalto that of the randomized reference system (Fig. 2(a)). In the RMS scenario, the overall shape of thep-value distribution is different from the BLS case, because RMS is a heterogeneous random walk.Nevertheless, the observed and reference distributions are again almost identical (Fig. 2(b)). In theTGS scenario, the shape of the p-value distribution is similar to that of RMS, because both searchstrategies share the feature of mode switching, one being controlled by chemoattractant gradients, theother occurring just randomly. Now, however, there are significant differences between the observedand reference distributions (Fig. 2(c)). In particular, the observed distribution shows a larger prob-ability of p-values smaller than 1/2, thus indicating attractive interactions. In the SGS scenario, wefind yet another shape of the p-value distribution, but again there are significant differences betweenthe observed and reference distributions (Fig. 2(d)). Application of p-value method to NK/K562 data
Next, we apply the p-value method to data from experiments with primary human NK immune cellsand K562 tumor cells in a three-dimensional collagen gel. We have evaluated nine independent datasets (see table 1 for details), using in each case a field of view of 867 µm × µm and a recordingtime interval of 0 .
75 min. Among the data sets, the total recording time T rec varied between 4 .
18 and19 .
35 hours, but in all cases a sufficiently large number N tri of triplets could be used for evaluation(See 4.th column of table 1. Note that N tri can exceed N fra · N (0) imm if new immune cells enter thefield of view over the recording period). The initial number N (0) imm of immune cells in the field of view7S N fra T rec (h) N tri N (0) imm N (0) tum v = σ ∆ t ( µmmin ) κ DS = 0 . . .
8. Here, N fra is the number of videoframes, T rec the total recording time in hours, N tri the number of valid triplets that could be used forthe p-value evaluation, N (0) imm the initial number of immune cells, N (0) tum the initial number of tumorcells, v the average speed of immune cells in µ m/min, and κ is the average persistence parameter ofimmune cells.varied between 13 and 47, the initial number N (0) tum of tumor cells between 14 and 26.In all data sets, we observed numerous encounters between the two cell types (see video V1 of theSupplemental Material, which shows the labeled cell positions of data set 0 during the whole recordingperiod). Nevertheless, the observed p-value distributions are located within the confidence intervalsof the reference system (Fig. 3(a-i)), indicating that the immune cell find their targets by chance.We attribute this surprising finding mainly to the unusually high migration speed of the NK immunecells (average migration speeds v varied between 4 .
12 and 6 . µm/min among the data sets), whichhelps them to explore a large search space within a relatively short time. Additionally, the NKimmune cells are able to migrate with high directional persistence (average persistence parameters κ varied between 1 .
73 and 3 . t = 0 and usedit as the initial configuration for a simulation of a blind random walk. Moreover, the parametersof the simulation were adjusted such that we obtained the same average migration properties ( v =5 . µm/min and κ = 2 .
59) as in the experiments. As a result, we indeed find that the simulatedblind migration produces a similar rate of immune-tumor cell encounters than seen in the data (seevideo V2 of the Supplemental Material).Finally, we performed another simulation, using again the initial configuration from data set 0 andmatching migration properties, but this time added the ability for spatial gradient sensing. Theresult demonstrates convincingly that, given their migration properties, chemotactical immune cellscould find the targets much faster than seen in the data (See video V3 of the Supplemental Material).
Effect of the recording interval ∆ t on p-value distributions Even if the motion of cells appears non-directional on short time scales, a target-directed migrationmay nevertheless emerge on larger time scales. To test for this possibility, we can sub-sample therecorded data by evaluating the triplets at time points t − n ∆ t , t , and t + n ∆ t , with an integer number n , thus effectively increasing the recording time period to n ∆ t . When applying this sub-samplingapproach to the surrogate data simulated in the SGS scenario (columns of Fig. 4), we indeed findthat the differences between the observed versus reference distributions become more pronounced8or larger effective recording intervals. By contrast, no strong effect is found after sub-sampling themeasured NK/K562 data (rows of Fig. 5). Generally, the reference distributions approach a uniformdistribution in the limit of very large recording intervals. Effect of the maximum detection radius r max on the p-value distributions We also test the effect of the maximum detection radius r max on the p-value distributions, both forthe SGS simulations (rows of Fig. 4) and for the measured data (columns of Fig. 5). Both resultsdemonstrate that a larger r max is generally preferable, because it reduces the widths of the confidenceintervals. When applying our method to new systems in which the range of interactions is unknown,we therefore recommend to set r max to the diagonal size of the field of view. Signature of weakly repulsive and weakly attractive interactions
Finally, we consider the case of very weak interactions between immune and target cells, using againthe simulated data in the SGS scenario. In order to modulate the interaction strength, we vary thechemotactic response parameter c (denoted by c R in [9]), which controls how sensitively the immunecells turn into the direction of the chemotactic concentration gradient. We find that for attractiveinteractions (positive c , bottom row of Fig. 6), p-values smaller than 1/2 are still more frequent thanin the reference system, but the differences eventually become non-significant in the case of veryweak attraction (case c = +5 in Fig. 6). Conversely, in the case of repulsive interactions (negative c ,top row of Fig. 6), p-values smaller than 1/2 are less frequent than in the reference system. Discussion
The collective motion of organisms can often be described in the framework of self-driven, interactingagents [16], provided the distance-dependent interactions between the agents are known. For thisreason, various methods have been developed in order to extract interaction functions from motiondata in groups of animals [17, 18, 19, 20, 21, 22, 23]To our knowledge, it has not yet been attempted to extract long-range interactions between differentcell types that presumably have a predator-prey relation. We have therefore tested different ap-proaches to detect and quantify remote interactions between immune and tumor cells, based solidlyon recorded cell trajectories.A first possible approach, used in some of the above studies [20, 21], is to set up an explicit model forthe migration and interaction of the agents, and then to fit the unknown model parameters directlyto the measured agent trajectories, for example using maximum likelihood optimization. We haveapplied such a method of parameter inference to simulated data of immune/tumor cell systems [24]and could correctly reproduce the known model parameters in some of the test cases. However, theinference produced wrong results whenever the investigated system had properties not fully capturedby the assumed migration and interaction model - unfortunately a common situation in biology.In this paper, we have therefore developed a new method of interaction detection based on p-values,which does not presume any detailed model of cell behavior, but only assumes that target-directedcells reveal themselves by a larger fraction of extraordinary, target-directed turns. Our methodhas only two user-adjusted parameters which slightly affect the results, namely the recording timeinterval ∆ t and the maximum expected interaction range r max . A third parameter, the minimumtrajectory length N ( tr ) min , does not usually need to be changed from its standard value 5. We havedemonstrated that this p-value method reliably distinguishes between target-blind migration andpurposeful pursuit in all test cases investigated so far.Recently, the misuse of p-values has been strongly criticized in the scientific community [25, 26, 27,28, 29]. The core of the problem is that many research studies treat the p-value as a uniquely definedfeature of their experiment, whereas there actually exists a (meta-) probability distribution for the9-value [30, 31]: When the very same experiment is repeated (that is, when new samples are drawnfrom the very same statistical model), the p-value will fall sometimes below and sometimes above thesignificance level. Picking just a single p-value thereby leads to non-reproducible results. For thisreason, we compute the probability distribution q obs (ˆ p ) of p-values (approximated by a histogramwith bins of finite size), pooled over all recorded steps of the immune cells. We then comparethis observed distribution with the distribution q ref (ˆ p ) of randomized reference systems. In orderto estimate the statistical fluctuation of the reference distribution, we compute q ref (ˆ p ) for a largeensemble of random reference systems, thus yielding a confidence interval for each bin of the p-valuehistogram. If there are interactions between immune and target cells, the observed p-values will befound outside of the confidence interval in at least some histogram bins.Our data sets with NK and K562 cells did not provide any evidence for long-range interactions,indicating that these cells meet each other only by chance. Although this particular negative resultdoes not rule out the existence of long-range interactions between other types of immune and tar-get cells, we speculate that the immune system might use an alternative way to effectively locatepathogens, which does not require any remote sensing abilities: Considering the huge variety ofpossible pathogens, and assuming that many of these pathogens do not advertise themselves to theimmune cells by emitting chemoattractant, the best method to increase the overall search efficiencyof the immune system is to make the immune cells large in number, fast and directionally persistent,thereby maximizing the rate of chance encounters with pathogens [13]. Combined with the abilitiesto recognize pathogens after steric contact and to recruit other immune cells by using endogenouschemoattractants, these properties create a robust and flexible system of defense.10 dditional information Author contributions statement:
CMe has devised the study, developed and implemented themethods for detecting interactions, applied the methods to the data, and wrote the paper. FH¨o hasperformed the measurements. CMa has developed and implemented the cell tracking algorithm. CBohas extracted, in vitro activated and expanded the immune cells. BFa has developed the imagingsystem and supervised the generation of the raw data.
Funding:
This work was funded by the Grant ME1260/11-1 (347962689) of the German ResearchFoundation DFG.
Competing interests statement:
The authors declare no competing interests.
Data availability statement:
All data and the Python implementation of the p-value method areavailable online at https://figshare.com/s/7d9a3dc191078289ba0e
Ethical approval and informed consent:
Not applicable.
Third party rights:
All material used in the paper are the intellectual property of the authors.11 eferences [1] Paul J Nahin.
Chases and escapes: the mathematics of pursuit and evasion . Princeton UniversityPress, 2012.[2] Howard C Berg.
Random walks in biology . Princeton University Press, 1993.[3] Michael Eisenbach.
Chemotaxis . World Scientific Publishing Company, 2004.[4] Charles A Janeway, Paul Travers, Mark Walport, Mark Shlomchik, et al.
Immunobiology: theimmune system in health and disease , volume 7. Current Biology London, 1996.[5] Himanshu Kumar, Taro Kawai, and Shizuo Akira. Pathogen recognition by the innate immunesystem.
International reviews of immunology , 30(1):16–34, 2011.[6] Manfred Schuster, Andreas Nechansky, and Ralf Kircheis. Cancer immunotherapy.
BiotechnologyJournal: Healthcare Nutrition Technology , 1(2):138–147, 2006.[7] Steven A Rosenberg. Decade in review—cancer immunotherapy: entering the mainstream ofcancer treatment.
Nature Reviews Clinical Oncology , 11(11):630, 2014.[8] Charles G Drake, Evan J Lipson, and Julie R Brahmer. Breathing new life into immunotherapy:review of melanoma, lung and kidney cancer.
Nature reviews Clinical oncology , 11(1):24, 2014.[9] Claus Metzner. On the efficiency of chemotactic pursuit - comparing blind search with temporaland spatial gradient sensing.
Scientific reports , 9(1):1–14, 2019.[10] Claus Metzner, Christoph Mark, Julian Steinwachs, Lena Lautscham, Franz Stadler, and BenFabry. Superstatistical analysis and modelling of heterogeneous random walks.
Nature commu-nications , 6(May):7516, jun 2015.[11] Frederick C Neidhardt. Escherichia coli and salmonella.
Typhimurium Cellular and MolecularBiology , 1987.[12] Xiao Zhou, Renping Zhao, Karsten Schwarz, Matthieu Mangeat, Eva C Schwarz, MohamedHamed, Ivan Bogeski, Volkhard Helms, Heiko Rieger, and Bin Qu. Bystander cells enhance nkcytotoxic efficiency by reducing search time.
Scientific reports , 7:44357, 2017.[13] Frederic Bartumeus, M G E da Luz, Gandhimohan M Viswanathan, and Jordi Catalan. Animalsearch strategies: a quantitative random-walk analysis.
Ecology , 86(11):3078–3087, 2005.[14] Merran Evans, Nicholas Hastings, and Brian Peacock. Statistical distributions. 2000.[15] Peter H Westfall, S Stanley Young, et al.
Resampling-based multiple testing: Examples andmethods for p-value adjustment , volume 279. John Wiley & Sons, 1993.[16] Tam´as Vicsek, Andr´as Czir´ok, Eshel Ben-Jacob, Inon Cohen, and Ofer Shochet. Novel type ofphase transition in a system of self-driven particles.
Physical review letters , 75(6):1226, 1995.[17] Michele Ballerini, Nicola Cabibbo, Raphael Candelier, Andrea Cavagna, Evaristo Cisbani, IreneGiardina, Vivien Lecomte, Alberto Orlandi, Giorgio Parisi, Andrea Procaccini, et al. Interactionruling animal collective behavior depends on topological rather than metric distance: Evidencefrom a field study.
Proceedings of the national academy of sciences , 105(4):1232–1237, 2008.[18] Ryan Lukeman, Yue-Xian Li, and Leah Edelstein-Keshet. Inferring individual rules from collec-tive behavior.
Proceedings of the National Academy of Sciences , 107(28):12576–12580, 2010.[19] Yael Katz, Kolbjørn Tunstrøm, Christos C Ioannou, Cristi´an Huepe, and Iain D Couzin. Infer-ring the structure and dynamics of interactions in schooling fish.
Proceedings of the NationalAcademy of Sciences , 108(46):18720–18725, 2011.1220] Anders Eriksson, Martin Nilsson Jacobi, Johan Nystr¨om, and Kolbjørn Tunstrøm. Determininginteraction rules in animal swarms.
Behavioral Ecology , 21(5):1106–1111, 2010.[21] Richard P Mann. Bayesian inference for identifying interaction rules in moving animal groups.
PloS one , 6(8):e22827, 2011.[22] Jacques Gautrais, Francesco Ginelli, Richard Fournier, St´ephane Blanco, Marc Soria, HuguesChat´e, and Guy Theraulaz. Deciphering interactions in moving animal groups.
Plos computa-tional biology , 8(9):e1002678, 2012.[23] Alessandro Attanasi, Andrea Cavagna, Lorenzo Del Castello, Irene Giardina, Tomas S Grigera,Asja Jeli´c, Stefania Melillo, Leonardo Parisi, Oliver Pohl, Edward Shen, et al. Informationtransfer and behavioural inertia in starling flocks.
Nature physics , 10(9):691, 2014.[24] Claus Metzner. Inferring long-range interactions between immune and tumor cells – pitfalls and(partial) solutions. arXiv preprint arXiv:1907.10284 , 2019.[25] Frank L Schmidt and John E Hunter. Eight common but false objections to the discontinuationof significance testing in the analysis of research data.
What if there were no significance tests ,pages 37–64, 1997.[26] Andrew Gelman and Hal Stern. The difference between significant and not significant is notitself statistically significant.
The American Statistician , 60(4):328–331, 2006.[27] Steven Goodman. A dirty dozen: twelve p-value misconceptions. In
Seminars in hematology ,volume 45, pages 135–140. Elsevier, 2008.[28] Valen E Johnson. Revised standards for statistical evidence.
Proceedings of the NationalAcademy of Sciences , 110(48):19313–19317, 2013.[29] Demetrios N Kyriacou. The enduring evolution of the p value.
Jama , 315(11):1113–1115, 2016.[30] Harold Sackrowitz and Ester Samuel-Cahn. P values as random variables - expected p values.
The American Statistician , 53(4):326–331, 1999.[31] Nassim Nicholas Taleb. A short note on p-value hacking. arXiv preprint arXiv:1603.07532 ,2016. 13 c)(b)(a) (1) (2) (3)
Figure 1: Explanation of the p-value based method to detect goal-directed migration. (a)
Compu-tation of the p-value (green shaded area under the curve p ( θ )). We consider a ’triplet’, consisting ofthree consecutive positions ( (cid:126)r ( i ) t − , (cid:126)r ( i ) t , (cid:126)r ( i ) t +1 ) of the focal immune cell (blue solid circles) and the posi-tion (cid:126)r ( i ) t of a target cell (red circle) in the vicinity. The ’persistence cone’ (right shaded area) is theinterval of the most probable migration directions of a target-blind immune cell, which is determinedby the previous migration direction (between time step t − t ), and by the known turning angledistribution p ( θ ) (olive curve on the top). The quantity θ (cid:63) (red) is the optimal turning angle thatwould align the immune cell precisely towards the target cell. The approach cone (left shaded area)is the interval of migration directions which are at least as target-oriented as the actual move of theimmune cell. By integrating p ( θ ) over the approach cone (that is, from turning angle θ = θ (cid:63) − ∆ θ to θ = θ (cid:63) +∆ θ ), we compute a p-value, the probability that a move at least as target-directed as observedcould occur in a target-blind random walk . (b) Three examples of immune cell trajectories (black)in relation to a target cell (red). Cases (1) and (2) are not indicative of goal-directed migration,but case (3) is ’suspicious’, because the cell steps out of its persistence cone and is at the same timevery target-directed. (c)
A reference system without interactions between immune and target cellsis generated by re-positioning the target cells randomly, while leaving the immune cell trajectoryunchanged. 14 a) (b)(c) (d)
Figure 2: Application of the p-value method to four different types of simulated data. Shown arein each case the p-value distributions q obs (ˆ p ) of the actual system (solid black lines) and the distri-butions q ref (ˆ p ) of the N s = 100 randomized reference systems (light gray lines). The colored dashedlines mark a confidence interval: values above the red dashed line (or below the blue dashed line)occur with a probability of 5 percent in the randomized system. (a) Blind search (BLS): Simulatedimmune cells migrate blindly, according to a correlated random walk with temporally constant mi-gration parameters. (b)
Random mode switching (RMS): Simulated immune cells migrate blindly,according to a correlated random walk with temporally fluctuating migration parameters. (c)
Tem-poral gradient sensing (TGS): Simulated immune cells use temporal gradients of a chemo-attractantto pursue the target cells. (d)
Spatial gradient sensing (SGS): Simulated immune cells use spatialgradients of a chemo-attractant to pursue the target cells. In cases (a) and (b), where there are nointeractions between simulated immune and target cells, the actual p-value distributions are insidethe confidence intervals of the randomized system. In cases (c) and (d), chemotactic interactionsbetween simulated immune and target cells lead to significant differences between the observed andreference distributions. 15 mm ~30µm (a) (b) (c)(d) (e) (f)(g) (h) (i)(j) (k) (l) Figure 3: Top: Application of the p-value method to nine measured data sets (a-i). In all cases(a-i), the confidence intervals are extremely narrow and the measured data are inside or very closeto these confidence intervals. This result is consistent with the null-hypothesis that immune cells are’blind’ with respect to the target cells. Bottom: Sketch of the experimental setup (j,k) and exampleimage with tracked cells (l). The K562 tumor cells are labeled with blue numbers, the NK cells withred numbers. 16igure 4: Top: Application of the p-value method to data simulated in the SGS model, usingdifferent recording intervals ∆ t (columns) and maximum interaction distances r max (rows). The dif-ferences between the simulated and reference data are more pronounced for larger recording intervals.Increasing the maximum interaction distance helps to reduce the width of the confidence interval inthe reference distributions. 17 t = 0.75 min � t = 1.50 min � t = 2.25 min r m a x = � m r m a x = � m r m a x = � m Figure 5: Top: Application of the p-value method to data set 0 (case (a) of Fig. 3), using differentrecording intervals (columns) and maximum interaction distances (rows). In all cases, the confidenceintervals are extremely narrow and the measured data are inside or very close to these confidenceintervals. 18 =-15 c=-10 c= -5c= +5 c=+10 c=+15
Figure 6: Top: Application of the p-value method to simulated data (modified SGS model), wheretumor cells are for the immune cells weakly repulsive (top row, negative values of the chemotaxisresponse coefficient c ), or weakly attractive (bottom row, positive values of the chemotaxis responsecoefficient c ). Note that in the standard SGS model (Fig. 2(d)), the coefficient is c = +500. Forrepulsive interactions, there are fewer small p-values and more large p-values than in the referencesystems. 19 upplemental information Video material
To compare the experimentally observed cell behavior with the models of blind search and spatialgradient sensing, we provide three videos ( https://figshare.com/s/7d9a3dc191078289ba0e ).The video
V1.flv shows the tracked cells of our data set 0. The NK immune cells are shown as redcircles, the K562 tumor cells as blue circles. All cells are labeled with unique numbers. Once a tumorcell is visited by an immune cell (within a distance smaller then 30 µm ), the tumor cell is consideredas ’found’ and is subsequently colored in gray.The video V2.flv shows a simulation that starts with the same initial configuration as in data set 0.The simulated immune cells also migrate with the same average speed and directional persistence as inthe experiment (Tumor cells are assumed to by completely stationary for simplicity). This simulationassumes a target-blind random walk (blind search BLS) and produces a rate of encounters betweenimmune and tumor cells comparable to that in data set 0.The video