Physical constraints on accuracy and persistence during breast cancer cell chemotaxis

TL;DR

This study combines theory, simulations, and experiments to explore the physical limits of chemotactic accuracy and persistence in breast cancer cells, revealing inherent constraints on their movement performance.

Contribution

It introduces a random walk model linking chemotactic accuracy and persistence, demonstrating their mutual constraints and inherent limitations.

Findings

Chemotactic accuracy increases with gradient strength.

Persistence and speed are not significantly affected by gradient strength.

Chemotactic accuracy and persistence are mutually constrained.

Abstract

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.

Tables

Table 1. Table 1: Table of parameters and values used in cellular Potts model (CPM) simulations. See text for more detailed reasoning behind values.

Parameter	Value	Reason
Total time $T$	$9$ h	Experiments
Recording interval $\Delta t$	$15$ min	Experiments
Background concentration $c$	$2.5$ nM	Experiments
Concentration gradient $g$	$5$ nM/mm	Experiments
Relaxed cell area $A_0$	$400\mu {\text{m}}^2$	Experiments
Lattice spacing $\ell$	$2$ $\mu$m	$\sim$$100$ sites per cell
Cell-environment contact energy $\alpha$	$2$ $\mu$m-1	$\alpha \sim {\ell }^{-1}$
Area deviation energy $\lambda$	0.01 $\mu$m-1	${\lambda }^{-1/2}\ll A_0$
Simulation time $\tau$	100 s
Total receptor number $N$	10,000	CI saturation
Bias strength $ϵ$	56	Calibrated via CI
Persistence strength $\eta$	107	Calibrated via DP
Polarization memory decay rate $r$	0.01 s-1	$r\sim {\tau }^{-1}$

Equations

$$\text{CI}\equiv ⟨\cos \theta ⟩,$$

$$\text{DP}\equiv ⟨\frac{|\text{displacement}|}{\text{distance}}⟩.$$

$$u=\alpha L+\lambda (A-A_0{)}^2,$$

$$\mathcal{P}=\{\begin{array}{ll}{e}^{-\left(\Delta u-w\right)}& \Delta u-w>0\\ 1& \Delta u-w\le 0,\end{array}$$

$$w=\Delta \vec{x}\cdot \vec{p},$$

$$\frac{\Delta \vec{p}}{\tau }=r(-\vec{p}+\eta \Delta {\hat{x}}_{\tau }+ϵ\vec{q}).$$

$$\vec{q}=⟨(n_i-\bar{n}){\hat{r}}_i⟩,$$

$$P(\theta |{\theta }^{\prime })=\underset{\mathrm{bias}}{\underbrace{b\cos \theta }}+\underset{\mathrm{persistence}}{\underbrace{\frac{e^{p\cos (\theta -{\theta }^{\prime })}}{2\pi I_0(p)}}},$$

$$⟨r^2⟩$$

$$+2\left(2z^2-1\right)\left(1-e^{-\tilde{M}}\right)$$

$$\displaystyle\qquad\qquad\quad\,\,\,+2z^{2}\left(1-e^{-\tilde{M}}\right)^{2}\bigg{]},$$

$$⟨x⟩$$

$$\mathrm{CI}$$

$$\mathrm{DP}$$

$$\frac{\partial c_i}{\partial t}=D_i\cdot \nabla c_i$$

$$⟨r^2⟩$$

$$⟨x{⟩}^2$$

$$⟨r^2⟩$$

$$⟨x{⟩}^2$$

$${\mathrm{CI}}^2$$

$${\mathrm{DP}}^2$$

$${\mathrm{CI}}^2=1-2ϵ\frac{1-{\mathrm{DP}}^2}{{\mathrm{DP}}^2}$$

Full text

AUTHOR SUMMARY

One of the most ubiquitous and important cell behaviors is chemotaxis: the ability to move in the direction of a chemical gradient. Due to its importance, key aspects of chemotaxis have been quantified for a variety of cells, including the accuracy, persistence, and speed of cell motion. However, whether these aspects are mutually constrained is poorly understood. Can a cell be accurate but not persistent, or vice versa? Here we use theory, simulations, and experiments on cancer cells to uncover mutual constraints on the properties of chemotaxis. Our results suggest that accuracy and persistence are mutually constrained.

INTRODUCTION

Chemotaxis plays a crucial role in many biological phenomena such as organism development, immune system targeting, and cancer progression iglesias2008navigating ; roussos2011chemotaxis ; kim2013cooperative ; varennes2016sense . Specifically, recent studies indicate that chemotaxis occurs during metastasis in many different types of cancer roussos2011chemotaxis ; friedl2011cancer ; witsch2010roles ; woodhouse1997general ; wang2004differential ; shields2007autologous . 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 kim2013cooperative and others have been devoted to the intracellular biochemical processes involved in cell motion petrie2009random . 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 shi2013interaction . 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 funamoto2001role ; mouneimne2006spatial ; van2007biased ; kay2008changing (Fig 1A); although occasionally it is defined in terms of the ratio of distances traveled nelson1975chemotaxis or number of motile cells iellem2001unique ; mayr2002vascular ; fiedler2005vegf in the presence vs. absence of the gradient. Directional persistence petrie2009random (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 mccutcheon1946chemotaxis , length ratio gorelik2014quantitative , or straightness index codling2008random ), although recent work has pointed out advantages of using the directional autocorrelation time gorelik2014quantitative ; dang2013inhibitory . 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 varennes2017emergent . 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.

RESULTS

Quantifying accuracy, persistence, and speed

We measure accuracy using the chemotactic index (CI) funamoto2001role ; mouneimne2006spatial ; van2007biased ; kay2008changing

[TABLE]

where $\theta$ 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 [math] and $1$; whereas in response to a repellent, CI usually falls between $-1$ and [math]. $\text{CI}=1$ represents perfectly accurate chemotaxis in which cell displacement is parallel to the gradient direction (Fig 1B, top two examples), and $\text{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 mccutcheon1946chemotaxis ; gorelik2014quantitative ; codling2008random (Fig 1A),

[TABLE]

Note that this ratio goes by several names mccutcheon1946chemotaxis ; gorelik2014quantitative ; codling2008random , 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 [math] and $1$, and once again intuitive sense can be made of either limit. If $\text{{{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 $\to 0$ in the limit of an infinitely long non-persistent trajectory.

An alternative measure of persistence is the directional autocorrelation time $\tau_{\rm AC}=\int_{0}^{\infty}dt^{\prime}\ \langle\cos(\theta_{t+t^{\prime}}-\theta_{t})\rangle$, where $t^{\prime}$ is the time difference between two points in a trajectory, and the average is taken over all starting times $t$ gorelik2014quantitative ; dang2013inhibitory . 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 Fig S1 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 $\Delta 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-$\beta$ 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 roussos2011chemotaxis ; friedl2011cancer ; witsch2010roles ; woodhouse1997general ; wang2004differential ; shields2007autologous . Among these soluble factors, transforming growth factor-$\beta$ (TGF-$\beta$) is a key environmental cue for the invasion process roussos2011chemotaxis ; luwor2015single ; giampieri2009localized ; ikushima2010tgfbeta ; kleuser200817 . Therefore, we use TGF-$\beta$ 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-$\beta$ 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-$\beta$ 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-$\beta$.

To verify that a graded TGF-$\beta$ profile is generated in the center channel, we utilize 10kDa FITC-dextran, whose hydrodynamic radius ($2.3$ nm) is similar to that of TGF-$\beta$ (approximately $2.4$ nm venturoli2005ficoll ). 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-$\beta$. See Materials and methods for details.

First, we perform a control experiment with no TGF-$\beta$ 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 $\theta$ is uniform, the distribution of $\cos\theta$ 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 E, 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 petrie2009random .

Next, we expose cells to a TGF-$\beta$ 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-$\beta$ 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-$\beta$ gradient strengths, $g=0$, $1$, $5$, and $50$ nM/mm, in three separate experiments each (Fig 4A-C; the trajectories for all experiments and $g$ values are shown in Fig S2). We see in Fig 4A that, consistent with Fig 3, the CI is zero for the control and increases with gradient strength $g$. 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 4B, consistent with Fig 3, the DP slightly decreases with the gradient strength although the decrease is roughly within error bars. Finally, we see in Fig 4C that the increase in the speed is small, achieving a statistically significant difference with the control only at the largest gradient strength, and that the trend is not monotonic.

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 varennes2016sense . 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$ $\mu$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 van2007biased and epithelial cells in epidermal growth factor gradients ellison2016cell , respectively varennes2016sense . This suggests that the response of MDA-MB-231 cells to TGF-$\beta$ 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-C), 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 $\mu$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) graner1992simulation ; swat2012multi , a lattice-based simulation that has been widely used to model cell migration szabo2010collective ; kabla2012collective ; varennes2016collective (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 $\alpha$ and has a basal area $A_{0}$ from which it can fluctuate at an energetic cost $\lambda$. This gives the energy function

[TABLE]

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 szabo2010collective . Specifically, for a lattice with $S$ total sites, one update step occurs in a fixed time $\tau$ 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

[TABLE]

where $\Delta u$ is the change in energy associated with the attempt. The bias term is defined as

[TABLE]

where $\Delta\vec{x}$ is the change in the cell’s center of mass caused by the attempt, and $\vec{p}$ 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 $\tau$ according to

[TABLE]

The first term in Eq. 6 represents exponential decay of $\vec{p}$ at a rate $r$. Thus, $r^{-1}$ characterizes the polarization vector’s memory timescale. The second term causes alignment of $\vec{p}$ with $\Delta\hat{x}_{\tau}$ according to a strength $\eta$, where $\Delta\hat{x}_{\tau}$ is a unit vector pointing in the direction of the displacement of the center of mass in the previous time step $\tau$. Thus, this term promotes persistence because it aligns $\vec{p}$ in the cell’s previous direction of motion. The third term causes alignment of $\vec{p}$ with $\vec{q}$ according to a strength $\epsilon$, where $\vec{q}$ contains the gradient sensing information, as defined below. Thus, this term promotes bias of motion in the gradient direction.

The sensing vector $\vec{q}$ is an abstract representation of the cell’s internal gradient sensing network and is defined as

[TABLE]

where the average is taken over all lattice sites $i$ that comprise the cell, and receptor saturation is incorporated as described below. The unit vector $\hat{r}_{i}$ points from the cell’s center of mass to site $i$, the integer $n_{i}$ represents the number of TGF-$\beta$ molecules detected by receptors at site $i$, and $\bar{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-$\beta$ 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-$\beta$ molecules in the vicinity of site $i$, which is sampled from a Poisson distribution whose mean is $(c+gx_{i})\ell^{3}$, where $\ell$ 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 $\vec{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, then $\vec{q}$ will point in that direction. Adaptive sensing has been observed in the TGF-$\beta$ pathway frick2017sensing in the form of fold-change detection shoval2010fold (for shallow gradients, subtraction as in Eq. 7 is similar to taking a ratio as in fold-change detection ellison2016cell ).

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 $\Delta 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, $\Delta t=15$ min, $c=2.5$ nM, and $g=5$ nM/mm are taken from the experiments. We estimate $A_{0}=400$ $\mu$m2 from the experiments, and we take $\ell=2$ $\mu$m, such that a cell typically comprises $A_{0}/\ell^{2}=100$ lattice sites. We find that realistic cell motion is sensitive to $\alpha$: when $\alpha$ is too small the cell is diffuse and unconnected, whereas when $\alpha$ is too large the cell does not move because the cost of perturbing the perimeter is too large. The crossover occurs around $\alpha\sim\ell^{-1}$ as expected, and therefore we set $\alpha$ on this order, to $\alpha=2$ $\mu$m*-1*. In contrast, we find that cell motion is not sensitive to $\lambda$ (apart from $\lambda=0$ for which the cell evaporates), and therefore we set $\lambda=0.01$ $\mu$m*-4* corresponding to typical area fluctuations of $\lambda^{-1/2}/A_{0}=2.5\%$. In order for our Poisson sampling procedure to be valid, the time step $\tau$ must be much larger than the timescale $\ell^{2}/D$ for an attractant molecule or receptor to diffuse with coefficient $D$ across a lattice site. Taking $D\sim 10$ $\mu$m2/s, we find $\tau\gg 0.4$ s. At the other end, we must have $\tau<\Delta t=900$ s for meaningful data collection. We find that within these bounds, results are not sensitive to $\tau$, and therefore we set $\tau$ on the larger end at $\tau=100$ s to reduce computational run time.

The parameters $N$, $\eta$, and $\epsilon$ are calibrated from the experimental data in Fig 4A-C. 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-$\beta$ receptors per cell wakefield1987distribution ; mitchell1992characterization , places the saturation level at roughly $g=50$ nM/mm as in the experiments (Fig 4D). We set $\epsilon=56$ $\mu$m*-1* and $\eta=107$ $\mu$m*-1* to calibrate their cognate observables, CI and DP, respectively, to the corresponding experimental values at $g=5$ nM/mm (Fig 4D and E).

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, $\tau$), 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}=\tau$. We conclude that the memory timescale of MDA-MB-231 cells is very short when responding to TGF-$\beta$ 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 F), i.e., tens of microns per hour. Second, we find that the speed shows little dependence 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 alt1980biased ; othmer1988models , 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 fraley2012dimensional . 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/\Delta t=36$ as in the experiments. The probability of a step making an angle $\theta$ with respect to the gradient direction is

[TABLE]

where $\theta^{\prime}$ 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 ($\theta=0$) and therefore promotes bias in that direction. It integrates to zero over its range ($-\pi<\theta<\pi$) 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 $-\pi<\theta<\pi$) whose sharpness grows with $p$. It is maximal at the previous angle $\theta^{\prime}$ and therefore promotes persistence. The normalization factor $I_{0}$ is the zeroth-order modified Bessel function of the first kind.

The requirement that $P(\theta|\theta^{\prime})$ be non-negative over the entire range of $\theta$ 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\pi)^{-1}$, and non-negativity requires $b<(2\pi)^{-1}\approx 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 othmer1988models . Specifically, the mean squared displacement and the mean displacement in the gradient direction are, in units of the mean step length,

[TABLE]

respectively, where $\tilde{M}=M(1-\psi)$ and $z=\chi/(1-\psi)$, with $\chi=\int_{-\pi}^{\pi}d\theta\ b\cos^{2}\theta=\pi b$ and $\psi=\int_{-\pi}^{\pi}d\phi\ [2\pi I_{0}(p)]^{-1}e^{p\cos\phi}\cos\phi=I_{1}(p)/I_{0}(p)$. We approximate the CI and DP in terms of these moments,

[TABLE]

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). In the limit $p=0$, Eqs. 11 and 12 are functions of only $b$ and $M$, and $b$ can be eliminated to yield ${\rm{{DP}}}=[1+M(1-{\rm CI}^{2})/2]^{-1/2}$ (left black line in Fig 6), where we have used the approximation $M\gg 1$ (see Materials and methods). Note here that when CI $=0$ we have DP $\approx(M/2)^{-1/2}$ for large $M$, which makes sense because for a simple random walk ($p=b=0$) the displacement goes like $M^{1/2}$ while the distance goes like $M$, such that DP $\sim M^{-1/2}$. Finally, the right edge corresponds to the maximal value of $p$ for a given $b$, for which we compute the approximation curve parametrically (right black line in Fig 6; see Materials and methods). We see in Fig 6 that these approximate expressions slightly underestimate the CI and overestimate the DP, but otherwise capture the crescent shape well. The under- and overestimation are due to the approximation $\langle r\rangle\approx\sqrt{\langle r^{2}\rangle}$ in Eqs. 11 and 12: because $\sigma_{r}^{2}=\langle r^{2}\rangle-\langle r\rangle^{2}\geq 0$ for any statistical quantity, we have $\sqrt{\langle r^{2}\rangle}\geq\langle r\rangle$, 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 $\to 1$, DP $\to 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 bias-induced 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-$\beta$ 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 dynamics, and more generally, it suggests that entire regions of the accuracy–persistence space are prohibited.

The present results provide some insights into TGF-$\beta$ induced migration mechanisms. Multiple signaling pathways induced by TGF-$\beta$ affect the dynamics of actin polymerization regulating cell migratory behaviorsikushima2010tgfbeta ; derynck2003smad-dependent ; moustakas2008dynamic ; olson2010linking . 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-$\beta$ olson2010linking ; dumont2003autocrine . 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 krause2014steering ; swaney2010eukaryotic ; edlund2002transforming ; fukata2003roles . The PI3K activity, thus, can possibly explain the present chemotactic responses of the breast cancer cells to TGF-$\beta$ 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 andrew2007chemotaxis ; bosgraaf2008pi3 . 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 krause2014steering ; swaney2010eukaryotic . The directional persistence could be more relevant to the polarity stability which is hardly controlled by chemotaxis krause2014steering as presented in the present results. In TGF-$\beta$ molecular cascades, activation of SMAD proteins could also affect the actin dynamics. Since SMAD-cascades include negative feedback inhibiting Rho activityderynck2003smad-dependent ; moustakas2008dynamic , 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 shi2013interaction . 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 shi2013interaction .

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 mugler2016limits ; ellison2016cell . 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.

MATERIALS AND METHODS

Cell culture and reagents

Human breast adenocarcinoma cells (MDA-MB-231) were cultured in Dulbecco’s Modified Eagle Medium/Ham’s F-12 (Advanced DMEM/F-12, Lifetechnologies, CA, USA) supplemented by 5$\%$ v/v fetal bovin serum (FBS), 2 mM L-glutamine (L-glu), and 100 $\mu$g ml-1 penicillin/streptomycin(P/S) for less than 15 passages. MDA-MB-231 cells were regularly harvested by 0.05$\%$ trypsin and 0.53mM EDTA (Lifetechnologies, CA, USA) when grown up to around 80$\%$ confluency in 75 cm2 T-flasks at 37 *∘*C with 5$\%$ CO2 incubation. Harvested cells were used for experiments or sub-cultured.

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$\%$ CO2. 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 $\mu$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 $\mu$g ml*-1* p/s rhee2010promigratory . After 24 hour-serum starvations, cells were exposed by a gradient of transforming growth factor beta-1 (TGF-$\beta$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 $\mu$m in thickness as described previously shin2013development . 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 $\mu$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:

[TABLE]

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 (FITC-dextran). 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 ozcelikkale2017effects ; ozcelikkale2017differential ; shin2016characterization , so that the microfluidic platform could be maintained at 37 *∘*C in a 5% CO2 environment during imaging. The time-lapse imaging was started 3 hours after applying TGF-$\beta$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 harley2008microarchitecture . Moreover, stationary cells due to the presence of the matrix were excluded harley2008microarchitecture ; giampieri2009localized ; clark2015modes ; haessler2012migration . 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-E 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 $\pm$ standard estimated error (S.E.M.) in Fig 4A-C. 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. Theoretical model reveals performance constraints and 10 become

[TABLE]

where $z=\pi b$, and we have neglected the exponential terms in the limit $M\gg 1$. Defining the small parameter $\epsilon=1/M$, these expressions become

[TABLE]

to first order in $\epsilon$, where $c\equiv 2(z^{-2}-2)$. Inserting these expressions into Eqs. 11 and 12, we obtain

[TABLE]

to first order in $\epsilon$. Because $z$ and $c$ are both functions only of $b$, we eliminate $b$ from Eqs. 18 and 19 to obtain

[TABLE]

to first order in $\epsilon$. 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\pi 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.

Acknowledgments

This work was supported by the Ralph W. and Grace M. Showalter Research Trust, the Purdue Research Foundation, and the Purdue University Center for Cancer Research Challenge Award.

Supporting information