Physical constraints on accuracy and persistence during breast cancer cell chemotaxis

Directed cell motion in response to an external chemical gradient occurs in many biological phenomena such as wound healing, angiogenesis, and cancer metastasis. Chemotaxis is often characterized by the accuracy, persistence, and speed of cell motion, but whether any of these quantities is physically constrained by the others is poorly understood. Using a combination of theory, simulations, and 3D chemotaxis assays on single metastatic breast cancer cells, we investigate the links among these different aspects of chemotactic performance. In particular, we observe in both experiments and simulations that the chemotactic accuracy, but not the persistence or speed, increases with the gradient strength. We use a random walk model to explain this result and to propose that cells’ chemotactic accuracy and persistence are mutually constrained. Our results suggest that key aspects of chemotactic performance are inherently limited regardless of how favorable the environmental conditions are.


Introduction
Chemotaxis plays a crucial role in many biological phenomena such as organism development, immune system targeting, and cancer progression [1][2][3][4]. Specifically, recent studies indicate that chemotaxis occurs during metastasis in many different types of cancer [2,[5][6][7][8][9]. At the onset of metastasis, tumor cells invade the surrounding extracellular environment, and oftentimes chemical signals in the environment can direct the migration of invading tumor cells. Several recent experiments have quantified chemotaxis of tumor cells in the presence of different chemoattractants [3] and others have been devoted to the intracellular biochemical processes involved in cell motion [10]. Since the largest cause of death in cancer patients is due to the metastasis, it is important to understand and prevent the directed and chemotactic behavior of invading tumor cells.
Chemotaxis requires sensing, polarization, and motility [11]. A cell's ability to execute these interrelated aspects of chemotaxis determines its performance. High chemotactic performance can be defined in terms of several properties. Cell motion should be accurate: cells should move in the actual gradient direction, not a different direction. Cell motion should be persistent: cells should not waste effort moving in random directions before ultimately drifting in the correct direction. Cell motion should be fast: cells should arrive at their destination in a timely manner.
Indeed, most studies of chemotaxis use one or more of these measures to quantify chemotactic performance. Accuracy is usually quantified by the so-called chemotactic index (CI), most often defined in terms of the angle made with the gradient direction [12][13][14][15] (Fig 1A); although occasionally it is defined in terms of the ratio of distances traveled [16] or number of motile cells [17][18][19] in the presence vs. absence of the gradient. Directional persistence [10] (DP) is usually quantified by the ratio of the magnitude of the cell's displacement (in any direction) to the total distance traveled by the cell (Fig 1A; sometimes called the McCutcheon index [20], length ratio [21], or straightness index [22]), although recent work has pointed out advantages of using the directional autocorrelation time [21,23]. Speed is usually quantified in terms of instantaneous speed along the trajectory or net speed over the entire assay.
However, the relationship among the accuracy, persistence, and speed in chemotaxis, and whether one quantity constrains the others, is not fully understood. Are there cells that are accurate but not very persistent, or persistent but not very accurate (Fig 1B)? If not, is it because such motion is possible but not fit, or is it because some aspect of cell motion fundamentally prohibits this combination of chemotactic properties?
Here we focus on how a cell's intrinsic migration mechanism as well as properties of the external environment place constraints on its chemotactic performance. The physics of diffusion places inherent limits on a cell's ability to sense chemical gradients [24]. These limits, along with the cell's internal information processing and its motility mechanism, determine the accuracy, persistence, and speed of migration. Using a human breast cancer cell line (MDA-MB-231) embedded within a 3D collagen matrix inside a microfluidic device imposing a chemical gradient, we are able to quantify the chemotactic performance of invasive cancer cells in response to various chemical concentration profiles. Results from chemotaxis assays are then compared with simulations and theoretical predictions in order to probe the physical limits of cancer cells to chemotaxis.

Quantifying accuracy, persistence, and speed
We measure accuracy using the chemotactic index (CI) [12][13][14][15] where θ is the angle the cell's displacement makes with the gradient direction (Fig 1A), and the average is taken over many cell trajectories. CI is bounded between −1 and 1. For chemotaxis in response to an attractant, as in this study, CI generally falls between 0 and 1; whereas in response to a repellent, CI usually falls between −1 and 0. CI = 1 represents perfectly accurate chemotaxis in which cell displacement is parallel to the gradient direction ( Fig 1B, top two examples), and CI = 0 indicates that the cells' migration is unbiased (Fig 1B, bottom two examples). The facts that CI is bounded and dimensionless make it easy to compare different values across different experimental conditions, and get an intuitive picture for the type of cell dynamics it represents. We measure persistence using the directional persistence (DP), defined as the ratio of the magnitude of the cell's displacement (in any direction) to the total distance traveled [20][21][22] (Fig 1A), Note that this ratio goes by several names [20][21][22], and although the name we use here contains the word 'chemotactic,' the ratio is in fact independent of the gradient direction. Indeed, DP measures the tendency of a cell to move in a straight line, in any direction. DP is also dimensionless and bounded between 0 and 1, and once again intuitive sense can be made of either limit. If DP = 1, then the cells are moving in perfectly straight lines in any arbitrary direction ( Fig 1B, right two examples). In contrast, a low DP is representative of a cell trajectory that starts and ends near the same location on average (Fig 1B, left two examples), with DP ! 0 in the limit of an infinitely long non-persistent trajectory.
An alternative measure of persistence is the directional autocorrelation time where t 0 is the time difference between two points in a trajectory, and the average is taken over all starting times t [21,23]. The advantage of the autocorrelation time is that, unlike the DP, it is largely independent of the measurement frequency and total observation time. The disadvantage is that, unlike the DP, it is not dimensionless or bounded. Although we use the DP here, we verify in S1 Fig that the autocorrelation time varies monotonically with the DP for our experimental assay.
We measure speed using the instantaneous speed along the trajectory. That is, we take the distance traveled in the measurement interval Δt (15 minutes in the experiments, see below), divide it by the interval, and average this quantity over all intervals that make up the trajectory.

Breast cancer cells chemotax up TGF-β gradients
We begin by investigating the above properties of chemotaxis in the context of metastasis, specifically the epithelial-mesenchymal transition and subsequent invasion of cancer cells. To this end, we perform experiments using a triple-negative human breast cancer cell line (MDA-MB-231). Invasion of tumor cells in vivo is aided by external cues including soluble factors that are thought to form gradients in the tumor microenvironment [2,[5][6][7][8][9]. Among these soluble factors, transforming growth factor-β (TGF-β) is a key environmental cue for the invasion process [2,[25][26][27][28]. Therefore, we use TGF-β as the chemoattractant.
The in vivo tumor microenvironment is highly complex. As a result, in vitro platforms have been developed and widely used to investigate the cancer response to a specific cue. In this study, a microfluidic platform is used to expose the TGF-β gradient to the cells in 3D culture condition (Fig 2A). The microfluidic device is designed with three different channels, a center, source, and sink channel ( Fig 2B). The center channel is filled with a composition of MDA-MB-231 cells and type I collagen while the medium is perfused through the side source and sink channels. TGF-β is applied only through the source channel, not the sink channel, and therefore a graded profile develops over time in the center channel by diffusion. Consequently, the MDA-MB-231 cells surrounded by type I collagen are exposed to a chemical gradient of TGF-β.
To verify that a graded TGF-β profile is generated in the center channel, we utilize 10kDa FITC-dextran, whose hydrodynamic radius (2.3 nm) is similar to that of TGF-β (approximately 2.4 nm [29]). The fluorescence intensity is shown in Fig 2C. The profile approaches steady state within 3 hours, is approximately linear, and remains roughly stationary for more than 12 hours. Therefore, we record the MDA-MB-231 cells using time-lapse microscopy every 15 minutes from 3 to 12 hours after imposing the TGF-β. See Materials and methods for details.
First, we perform a control experiment with no TGF-β to characterize the baseline of the MDA-MB-231 cell migratory behavior. Representative trajectories are shown in Fig 3A, and we see that there is no apparent preferred direction. Indeed, as seen in Fig 3C (black), the CI is centered around zero, indicating no directional bias. Notably, the spread of the CI values is very broad, with many data points falling near the endpoints −1 and 1. This is a generic feature of the CI due to its definition as a cosine: when the distribution of angles θ is uniform, the distribution of cos θ is skewed toward −1 and 1 because of the cosine's nonlinear shape. Nonetheless, we see that the median of the CI is very near zero as expected. The speed and DP are shown in Fig 3D and 3E, respectively (black). We see that the DP is significantly above zero, indicating that even in the absence of any chemoattractant, cells exhibit persistent motion. This result is consistent with previous works that showed that cells cultured in 3D tend to have directionally persistent movement unlike those in 2D [10].
Next, we expose cells to a TGF-β gradient of g = 50 nM/mm. Representative trajectories are shown in Fig 3B, and we see a possible bias in the gradient direction. Indeed, as seen in Fig 3C  (red), the CI is centered above zero, indicating a directional bias, and the difference with the control distribution is statistically significant (p value < 0.05). We also see in Fig 3D (red) that the speed increases, although we will see below that the increase is relatively small and that the trend is non necessarily monotonic. Finally, we see in Fig 3E (red) that the DP decreases, although the difference with the control is not statistically significant. These results suggest that a TGF-β gradient causes a significant increase in directional bias (CI) but not necessarily a significant change in cell speed or persistence (DP).
To confirm the trends suggested above, we evaluate the response to four different TGF-β gradient strengths, g = 0, 1, 5, and 50 nM/mm, in three separate experiments each ( In fact, the CI appears to saturate beyond 5 nM/mm, such that its value at 50 nM/mm is not significantly larger than its value at 5 nM/mm. We also see in Fig

Minimum detectable gradient is shallow
A striking feature of Fig 4A is that the cells respond to a gradient as shallow as g = 5 nM/mm. To put this value in perspective, we estimate both the relative concentration change and the absolute molecule number difference across the cell body [4]. The microfluidic device is about 1 mm in the gradient direction, and therefore a cell in the middle experiences a background concentration of about c = 2.5 nM. Assuming the cell is on the order of a = 10 μm wide, the change in concentration across its body is ga = 0.05 nM, for a relative change of ga/c = 2%. The number of attractant molecules that would occupy half the cell body is on the order of ca 3 = 1500. Two percent of this is ga 4 = 30, meaning that cells experience about a thirty-molecule difference between their two halves. The same quantities are approximately ga/c = 1% and 6%, and ga 4 = 60 and 300, for amoebae in cyclic adenosine monophosphate gradients [14] and epithelial cells in epidermal growth factor gradients [30], respectively [4]. This suggests that the response of MDA-MB-231 cells to TGF-β gradients is close to the physical detection limit for single cells.

Simulations suggest sensing and persistence are decoupled
To understand the experimental observation that the CI increases with gradient strength, but the DP and speed do not (Fig 4A-4C), we turn to computer simulations. The cells in the experiments are executing 3D migration through the collagen matrix (as opposed to crawling on top of a 2D substrate). Nevertheless, the imaging is acquired as a 2D projection of the 3D motion. We do not expect this projection to introduce much error into the analysis because the height of the microfluidic device is less than 100 μm, whereas its width in the gradient direction is about 1 mm, and its length is several millimeters. Indeed, from the experimental trajectories (Fig 3) we have estimated that if motility fluctuations in the height direction are equivalent to those in the length direction, then the error in the CI that we make by the fact that we only observe a 2D projection of cell motion is less than 1%. Consequently, for simplicity we use a 2D rather than 3D simulation of chemotaxis of a cell through an extracellular medium.
Specifically, we use the cellular Potts model (CPM) [31,32], a lattice-based simulation that has been widely used to model cell migration [33][34][35] (note that whereas often the CPM is used to model collective migration, here we use it for single-cell migration). In the CPM, a cell is defined as a finite set of simply connected sites on a regular square lattice (Fig 5). The cell adheres to the surrounding collagen with an adhesion energy α and has a basal area A 0 from which it can fluctuate at an energetic cost λ. This gives the energy function where L and A are the cell's perimeter and area, respectively. Cell motion is a consequence of minimizing the energy u subject to thermal noise and a bias term w that incorporates the response to the gradient [33]. Specifically, for a lattice with S total sites, one update step occurs in a fixed time τ and consists of S attempts to copy a random site's label (cell or non-cell) to a randomly chosen neighboring site. Each attempt is accepted with probability where Δu is the change in energy associated with the attempt. The bias term is defined as where Dx is the change in the cell's center of mass caused by the attempt, andp is the cell's polarization vector (Fig 5 inset, black arrow), described below. The dot product acts to bias cell motion because movement parallel to the polarization vector results in a more positive w, and thus a higher acceptance probability (Eq 4). The polarization vector is updated every time step τ according to The first term in Eq 6 represents exponential decay ofp at a rate r. Thus, r −1 characterizes the polarization vector's memory timescale. The second term causes alignment ofp with Dx t according to a strength η, where Dx t is a unit vector pointing in the direction of the displacement of the center of mass in the previous time step τ. Thus, this term promotes persistence because it alignsp in the cell's previous direction of motion. The third term causes alignment ofp withq according to a strength �, whereq contains the gradient sensing information, as defined below. Thus, this term promotes bias of motion in the gradient direction. The sensing vectorq is an abstract representation of the cell's internal gradient sensing network and is defined asq where the average is taken over all lattice sites i that comprise the cell, and receptor saturation is incorporated as described below. The unit vectorr i points from the cell's center of mass to site i, the integer n i represents the number of TGF-β molecules detected by receptors at site i, and � n is the average of n i over all sites. The integer n i is the minimum of two quantities: (i) the number of TGF-β receptors at site i, which is sampled from a Poisson distribution whose mean is the total receptor number N divided by the number of sites; and (ii) the number of TGF-β molecules in the vicinity of site i, which is sampled from a Poisson distribution whose mean is (c + gx i )ℓ 3 , where ℓ is the lattice spacing, and x i is the position of site i along the gradient direction. Taking the minimum incorporates receptor saturation, since each site cannot detect more attractant molecules than its number of receptors. The subtraction in Eq 7 makes q a representation of adaptive gradient sensing: if receptors on one side of the cell detect molecule numbers that are higher than those on the other side, thenq will point in that direction. Adaptive sensing has been observed in the TGF-β pathway [36] in the form of fold-change detection [37] (for shallow gradients, subtraction as in Eq 7 is similar to taking a ratio as in fold-change detection [30]).
The simulation is performed at a fixed background concentration c and gradient g for a total time T. The position of the cell's center of mass is recorded at time intervals Δt, from which we compute the CI, DP, and speed.
The parameter values used in the simulation are listed in Table 1 and are set in the following way. The values T = 9 h, Δt = 15 min, c = 2.5 nM, and g = 5 nM/mm are taken from the experiments. We estimate A 0 = 400 μm 2 from the experiments, and we take ℓ = 2 μm, such that a cell typically comprises A 0 /ℓ 2 = 100 lattice sites. We find that realistic cell motion is sensitive to α: when α is too small the cell is diffuse and unconnected, whereas when α is too large the cell does not move because the cost of perturbing the perimeter is too large. The crossover occurs around α � ℓ −1 as expected, and therefore we set α on this order, to α = 2 μm −1 . In contrast, we find that cell motion is not sensitive to λ (apart from λ = 0 for which the cell evaporates), and therefore we set λ = 0.01 μm −4 corresponding to typical area fluctuations of λ −1/2 /A 0 = 2.5%. In order for our Poisson sampling procedure to be valid, the time step τ must be much larger than the timescale ℓ 2 /D for an attractant molecule or receptor to diffuse with coefficient D across a lattice site. Taking D � 10 μm 2 /s, we find τ � 0.4 s. At the other end, we must have τ < Δt = 900 s for meaningful data collection. We find that within these bounds, results are not sensitive to τ, and therefore we set τ on the larger end at τ = 100 s to reduce computational run time.
The parameters N, η, and � are calibrated from the experimental data in Fig 4A-4C. Specifically, N sets the gradient value above which the CI saturates (see Fig 4A) because if the gradient is large but N is small, the cell quickly migrates into a region in which there are more attractant molecules than receptors at all lattice sites, and gradient detection is not possible. We find that N = 10,000, which is a reasonable value for the number of TGF-β receptors per cell [38,39], places the saturation level at roughly g = 50 nM/mm as in the experiments (Fig 4D). We set � = 56 μm −1 and η = 107 μm −1 to calibrate their cognate observables, CI and DP, respectively, to the corresponding experimental values at g = 5 nM/mm (Fig 4D and 4E).
The final parameter is the memory timescale of the polarization vector, r −1 . As seen in Fig  4E (gray), we find that the behavior of the DP depends sensitively on this timescale. When r −1 is large, the DP increases with gradient strength. In contrast, when r −1 is small (indeed, equal to the smallest timescale in the system, τ), the DP does not increase with gradient strength, and in fact slightly decreases (Fig 4E, blue). Because the latter behavior is consistent with the experiments (Fig 4B), we set r −1 = τ. We conclude that the memory timescale of MDA-MB-231 cells is very short when responding to TGF-β gradients.
We validate the simulation in two ways, using the speed. First, we find that the magnitude of the speed in the simulations is on the same order as the speed in the experiments (Fig 4C  and 4F), i.e., tens of microns per hour. Second, we find that the speed shows little dependence Physical constraints of cancer cell chemotaxis on the gradient strength in both the simulations and the experiments: it slightly increases in Fig 4C and slightly decreases in Fig 4F. Considering that the speed is not calibrated directly in our simulations, these consistencies validate the CPM as a reasonable description of the cell migration in the experiments. Our finding that the cell's memory timescale r −1 takes its minimum value allows for the following interpretation: the parameter r couples the persistence term and the sensory term in the CPM (Eq 6). Thus, when the memory timescale r −1 is long, biased motion must be also persistent and vice versa. In contrast, when the memory timescale r −1 is short, it is possible for bias to increase without increasing persistence. Therefore, the simulations suggest that the reason that CI but not DP increases with gradient strength in the experiments, is that the drivers of sensory bias and migratory persistence in the cell's internal network are decoupled from one another.

Theoretical model reveals performance constraints
Our finding that bias and persistence are decoupled in the simulations allows us to appeal to a much more simplified theoretical model in order to understand and predict global constraints on chemotaxis performance. Specifically, we consider the biased persistence random walk (BPRW) model [40,41], in which bias and persistence enter as explicitly independent terms controlled by separate parameters. The BPRW has been shown to be sufficient to capture random and directional, but not periodic, behaviors of 3D cell migration [42]. Because we do not observe periodic back-and-forth motion of cells in our experiments, we propose that the BPRW is sufficient to investigate chemotactic constraints here.
As in the simulations, we consider the BPRW model in 2D. In the BPRW model, a cell is idealized as a single point. Its trajectory consists of M steps whose lengths are drawn from an exponential distribution. We take M = T/Δt = 36 as in the experiments. The probability of a step making an angle θ with respect to the gradient direction is where θ 0 is the angle corresponding to the previous step. The first term incorporates the bias, with strength b. It is maximal when the step points in the gradient direction (θ = 0) and therefore promotes bias in that direction. It integrates to zero over its range (−π < θ < π) because the bias term only reshapes the distribution without adding or subtracting net probability. The second term incorporates the persistence, with strength p. It is a von Mises distribution (similar to a Gaussian distribution, but normalized over the finite range −π < θ < π) whose sharpness grows with p. It is maximal at the previous angle θ 0 and therefore promotes persistence. The normalization factor I 0 is the zeroth-order modified Bessel function of the first kind. The requirement that P(θ|θ 0 ) be non-negative over the entire range of θ mutually constrains b and p. However, apart from this constraint, b and p can take any positive value. We sample many pairs of b and p, reject those that violate the constraint, and compute the CI and DP from a trajectory generated by each remaining pair. The results are shown in Fig 6 (colored circles). We see in Fig 6 that the BPRW model exists in a highly restricted 'crescent' shape within CI-DP space. As expected, the CI increases with the bias parameter b (color of circles, from blue to red). The top corner corresponds to maximal bias and no persistence; indeed, when p = 0 the persistence term in Eq 8 reduces to (2π) −1 , and non-negativity requires b < (2π) −1 � 0.16, which is consistent with the upper limit of the color bar. Also as expected, the DP increases with the persistence parameter p (size of circles, from small to large), although only in the lower portion where the CI is low.
The crescent shape of the allowed CI and DP values in Fig 6 can be understood quantitatively because several moments of the BPRW are known analytically [41]. Specifically, the mean squared displacement and the mean displacement in the gradient direction are, in units of the mean step length, respectively, whereM ¼ Mð1 À cÞ and z = χ/(1 − ψ), with w ¼ R p À p dy b cos 2 y ¼ pb and c ¼ R p À p d� ½2pI 0 ðpÞ� À 1 e p cos � cos � ¼ I 1 ðpÞ=I 0 ðpÞ. We approximate the CI and DP in terms of these moments, and evaluate these expressions in specific limits to approximate the edges of the shape. In the limit b = 0, Eq 11 reduces to CI = 0 (bottom black line in Fig 6). because s 2 r ¼ hr 2 i À hri 2 � 0 for any statistical quantity, we have ffi ffi ffi ffi ffi ffi ffi ffi hr 2 i p � hri, making Eq 11 an underestimate and Eq 12 an overestimate.
The crescent shape can also be understood intuitively. First, we see that the DP cannot be smaller than a minimum value (region I in Fig 6). This is because the trajectory length M is finite, and as discussed above, the DP only vanishes for infinitely long trajectories. If M were to increase, the crescent would extend further toward DP = 0. Second, we see that the top of the crescent bends away from the CI !1, DP !0 corner (region II in Fig 6). In other words, it is not possible to have high bias without any persistence. This is because if the bias is strong, then cells will track the gradient very well. Consequently, they will move in nearly straight lines in the gradient direction, and straight movement corresponds to high persistence. This is a biasinduced persistence, distinct from the bias-independent persistence in the lower-right corner of the crescent. Finally, we see that the bending shape of the crescent implies that no solutions exist at large DP and intermediate CI (region III in Fig 6). In other words, it is not possible to have high persistence with partial bias. This is because, as mentioned above, persistence is induced either (i) directly, as a result of a large persistence parameter p which is independent of the bias, in which case the CI is low; or (ii) indirectly, as a result of a large bias parameter b, in which case the CI is high. Neither of these mechanisms permits intermediate bias, and therefore high persistence can be accompanied only by low or high directionality. Together, these features of the crescent shape imply that specific modes of chemotaxis are prohibited under our simple model, as indicated by the regions I, II, and III.
Finally, the crescent shape provides a qualitative rationale for the data from the simulations and experiments, which are overlaid in the cyan and red squares in Fig 6, respectively. Specifically, the shape of the crescent is such that if a cell has a low CI and intermediate DP (bottom right corner of the crescent) and its CI increases, its DP must decrease (solid magenta arrow in Fig 6). In contrast, a simultaneous increase in CI and DP from this starting position is not possible according to the model (dashed magenta arrow in Fig 6). We see that the data are qualitatively consistent with this predicted trend, as an increase in the CI corresponds to a decrease in the DP in both the experiments and the simulations (Fig 6, squares). There is quantitative disagreement, in the sense that the data do not quite overlap with the crescent, but this is a reflection of the extreme simplicity of the BPRW model. Nonetheless, the qualitative features of the BPRW model are sufficient to explain the way in which accuracy and persistence are mutually constrained during the chemotaxis response of these cells.

Discussion
By integrating experiments with theory and simulations, we have investigated mutual constraints on the accuracy (CI), persistence (DP), and speed of cancer cell motion in response to a chemical attractant. We have found that while the CI of breast cancer cells increases with the strength of a TGF-β gradient, the speed does not show a strong trend, and the DP slightly decreases. The simulations suggest that the decrease in DP is due to a decoupling between sensing and persistence in the migration dynamics. The theory confirms that the decrease in DP is due to a mutual constraint on accuracy and persistence for this type of decoupled Physical constraints of cancer cell chemotaxis dynamics, and more generally, it suggests that entire regions of the accuracy-persistence space are prohibited.
The present results provide some insights into TGF-β induced migration mechanisms. Multiple signaling pathways induced by TGF-β affect the dynamics of actin polymerization regulating cell migratory behaviors [27,[43][44][45]. Among these, phosphatidylinositol 3-kinase (PI3K) and the small GTPase-Rac1 signaling have been reported to promote actin organization of breast cancer cells in response to TGF-β [45,46]. PI3K and the Rho-family GTPase networks (including Rac1, RhoA and Cdc42) have been widely studied in chemotaxis, which regulates cell polarity and directional sensing [47][48][49][50]. The PI3K activity, thus, can possibly explain the present chemotactic responses of the breast cancer cells to TGF-β gradient. Recent studies have shown that PI3K is relevant to the accuracy of the cell movement in shallow chemoattractants, whereas it does not induce the orientation of cell movement in steep gradients; rather, PI3K contributes the motility enhancement [51,52]. These results can be correlated with the cell motility trend in the present experimental results. In addition, the PI3K signaling pathway has been reported not to mediate the persistence of cell protrusions which could be directly related to the DP [47,48]. The directional persistence could be more relevant to the polarity stability which is hardly controlled by chemotaxis [47] as presented in the present results. In TGF-β molecular cascades, activation of SMAD proteins could also affect the actin dynamics. Since SMAD-cascades include negative feedback inhibiting Rho activity [43,44], it may affect the cell responses highly promoted in CI but not in speed. However, the underlying molecular mechanisms need further research.
Our finding that sensing and persistence are largely decoupled in the migration dynamics is related to the view that directional sensing and polarity are separate but connected modules in chemotaxis [11]. Indeed, CI, DP, and speed in our study play the roles of the directional sensing, polarity, and motility modules, respectively, that have been shown to reproduce many of the observed behaviors of chemotaxing cells. Moreover, several of the the molecular signaling pathways discussed above, including those involving PI3K and Rho family GTPases, have been proposed as the potential networks corresponding to these modules [11].
Several predictions arise from our work that would be interesting to test in future experiments. First, our simulation scheme assumes that the saturation of the CI with gradient strength (Fig 4A) is due to limited receptor numbers. However, alternative explanations exist that are independent of the receptors, such as the fact that it is more difficult to detect a concentration difference on top of a large concentration background than on top of a small concentration background due to intrinsic fluctuations in molecule number [30,53]. An interesting consequence of our mechanism of receptor saturation is that, at very large gradients (beyond those of Fig 4A), the CI would actually decrease because all receptors would be bound. It would be interesting to test this prediction in future experiments.
Second, our work suggests that not all quadrants of the accuracy-persistence plane are possible for cells to achieve (Fig 6). It would be interesting to measure the CI and DP of other cell types, in other chemical or mechanical environments, to see if the crescent shape seen in Fig 6  is a universal restriction, or if not, what new features of chemotaxis are therefore not captured by the modeling. In this respect, the work here can be seen as a null model, deviations from which would indicate new and unique types of cell motion.
Cell-matrix composition was prepared in the microfluidic device. For the composition, MDA-MB-231 cells were mixed with 2 mg/ml of type I collagen (Corning Inc., NY, USA) mixture prepared with 10X PBS, NaOH, HEPE solution, FBS, Glu, P/S, and cell-culture level distilled water after centrifuged with 1000 rpm for 3 minutes. The cell mixture was filled in center-channel of the microfluidic devices and incubated in at 37˚C with 5% CO 2 . The cells in the collagen matrix were initially cultured in basic medium (DMEM/F12 supplemented by 5% v/v FBS, 2 mM L-glu, and 100 μg ml −1 p/s) for 24 hours. Then the cells were exposed by reduced serum medium for another 24 hours, which was advanced DMEM/F12 containing 1% v/v FBS, 2 mM L-glu, and 100 μg ml −1 p/s [54]. After 24 hour-serum starvations, cells were exposed by a gradient of transforming growth factor beta-1 (TGF-β1, Invitrogen, CA, USA).

Microfluidic device for chemical gradient
The microfluidic device was designed to generate a linear gradient of soluble factors (Fig 2). The device is composed of three channels which are 100 μm in thickness as described previously [55]. A center channel that is 1 mm wide aims to culture tumor cells with ECM components. The center channel is connected to two side channels. The 300 μm-wide side channels are connected to large reservoirs at the end ports including culture medium. Since the side channels are in contact with the top and bottom sides of the center channel, the growth factor gradient can be generated by diffusing the soluble factor from one of the side channels, a source channel, to the other, a sink channel. Assuming there is neither pressure difference nor flow between the side channels, the concentration of a given factor can be described by the chemical species conservation equation as follows: Once the concentration profile in the center channel reaches steady state, the linear profile persists for a while and can therefore be approximated by assuming the boundary conditions of concentration at the side channels are constants. To verify the diffusion behavior, the gradient formation was examined by using 10k Da FITC-fluorescence conjugated dextran (FITCdextran). FITC-dextran solution was applied in the source channel while the sink channel was filled with normal culture medium. The FITC-dextran concentration profile was evaluated by the FITC fluorescent intensity in the center channel. To disregard the effect of photo-bleaching on the results, the intensity was normalized by the intensity of the source channel. The normalized intensity was reasonably considered since the fluorescence intensity of the source channel consistently remained as maximum due to the large reservoirs. The FITC dextran intensity profile (Fig 2C) showed that the linear profile was developed within 3 hours after applying the source and continued for more than 9 hours.

Characterization of cell migration with time-lapse microscopy
Cell behaviors were captured every 15 minutes for 9 hours using an inverted microscope (Olympus IX71, Japan) equipped with a stage top incubator as described previously [56][57][58], so that the microfluidic platform could be maintained at 37˚C in a 5% CO 2 environment during imaging. The time-lapse imaging was started 3 hours after applying TGF-β1 solution in the source channel to have sufficient adjusting time. To analyze each cell behavior, a cell area in the bright field images were defined by a contrast difference between the cells and a background, and the images were converted to monochrome images by using ImageJ. Cell trajectories were demonstrated by tracking centroids of the cell area. In tracking the cell movements, cells undergoing division were excluded to avoid extra influences to affect cell polarity [59]. Moreover, stationary cells due to the presence of the matrix were excluded [26,[59][60][61]. The stationary cells were defined as the cells that moved less than their diameter. A migration trajectory was defined by connecting the centroids of a cell from each time point.

Statistical analysis of experiments
In examining the chemotactic characteristics of each group, more than 40 cell trajectories were evaluated per a group. A data point in Fig 3C-3E indicates each metric of a cell trajectory showing distribution characteristics with a box plot. The box plot includes boundaries as quadrants and a center as a median. The distribution of each metric was statistically analyzed by using Mann-Whitney U-test. This non-parametric method was used since the distribution was not consistently normal (the CI is a function of cosine). The significant change on the population lies on the biased distribution of each cell parameter when the p value < 0.05. Furthermore, the experiments were repeated at least 3 times and reported with means of medians ± standard estimated error (S.E.M.) in Fig 4A-4C. To evaluate physical limits on each metric, the data points were compared each other using a student t-test. The statistical significance between comparisons were examined when the p value < 0.05.

Mathematical approximations
In the limit p = 0, Eqs 9 and 10 become hr 2 i ¼ z 2 M 2 þ 2ð1 À 2z 2 ÞM þ 2ð3z 2 À 1Þ; ð14Þ where z = πb, and we have neglected the exponential terms in the limit M � 1. Defining the small parameter � = 1/M, these expressions become to first order in �, where c � 2(z −2 − 2). Inserting these expressions into Eqs 11 and 12, we obtain to first order in �. Because z and c are both functions only of b, we eliminate b from Eqs 18 and 19 to obtain to first order in �. This expression is equivalent to that given below Eq 12 and provides the left black line in Fig 6. The right black line in Fig 6 corresponds to the maximal value of p for a given b that keeps Eq 8 non-negative. Non-negativity requires that the sum of the minimal values of each term in Eq 8 is zero: −b + e −p /[2πI 0 (p)] = 0. With this expression for b in terms of p, Eqs 11 and 12 become functions of only p and M. Therefore, by varying p, we compute the right black line parametrically.