Quantifying accuracy, persistence, and speed

We measure accuracy using the chemotactic index (CI) [12–15] (1) 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–22] (Fig 1A), (2) Note that this ratio goes by several names [20–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′ 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–9]. Among these soluble factors, transforming growth factor-β (TGF-β) is a key environmental cue for the invasion process [2, 25–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-β.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Microfluidic device used as a chemotaxis platform.

(A) Cross-sectional view illustrating concentration gradient formed by diffusion. (B) Illustration showing structure of the microfluidic channels. Center channel (green) is filled with type I collagen mixture and MDA-MB-231 mixture, source channel is filled with culture medium containing TGF-β, and sink channel is filled with only culture medium. (C) FITC-dextran fluorescence within the center channel. Blue region indicates sink channel while red region indicates source channel.

https://doi.org/10.1371/journal.pcbi.1006961.g002

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].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Cell trajectories and chemotaxis metrics.

Cell trajectories of (A) control and (B) 50nM/mm TGF-β gradient. Distribution of (C) chemotactic index, (D) speed, and (E) directional persistence of each trajectory from both the control (black) and the TGF-β gradient (red). Boundary of box plots indicates quadrants with centerline as median. Distributions are statistically compared using Mann-Whitney test.

https://doi.org/10.1371/journal.pcbi.1006961.g003

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 (Fig 4A–4C; the trajectories for all experiments and g values are shown in S2 Fig). 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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Comparison of experiments with simulations.

Experimental (A) chemotactic index, (B) directional persistence, and (C) speed for four different TGF-β gradients, g = 0, 1, 5, and 50 nM/mm.(red) Data points indicate average and standard error of medians from three different experiments. A, B, and C are plotted with log-scaled TGF-β gradient. (D-F) Same for cellular Potts model (CPM) simulations (blue). Error bars are standard error from 1000 trials. Directional persistence from reduced polarization memory decay rate(r) is represented in (E) (gray).

https://doi.org/10.1371/journal.pcbi.1006961.g004

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³ = 1500. Two percent of this is ga⁴ = 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⁴ = 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–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₀ from which it can fluctuate at an energetic cost λ. This gives the energy function (3) where L and A are the cell’s perimeter and area, respectively.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Cellular Potts model (CPM) simulation.

Snapshot shows cell (gray) migrating towards increasing chemical concentration over time (white trajectory). Inset: Cell motility occurs through addition and removal of lattice sites. , cell polarization vector; α, cell-collagen adhesion energy; δL, change in perimeter; δA, change in area.

https://doi.org/10.1371/journal.pcbi.1006961.g005

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 (4) where Δu is the change in energy associated with the attempt. The bias term is defined as (5) where is the change in the cell’s center of mass caused by the attempt, and 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 (6) The first term in Eq 6 represents exponential decay of at a rate r. Thus, r⁻¹ characterizes the polarization vector’s memory timescale. The second term causes alignment of with according to a strength η, where 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 aligns in the cell’s previous direction of motion. The third term causes alignment of with according to a strength ϵ, where contains the gradient sensing information, as defined below. Thus, this term promotes bias of motion in the gradient direction.

The sensing vector is an abstract representation of the cell’s internal gradient sensing network and is defined as (7) 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 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 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)ℓ³, 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 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 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₀ = 400 μm² from the experiments, and we take ℓ = 2 μm, such that a cell typically comprises A₀/ℓ² = 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 α ∼ ℓ⁻¹ as expected, and therefore we set α on this order, to α = 2 μm⁻¹. 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⁻⁴ corresponding to typical area fluctuations of λ^−1/2/A₀ = 2.5%. In order for our Poisson sampling procedure to be valid, the time step τ must be much larger than the timescale ℓ²/D for an attractant molecule or receptor to diffuse with coefficient D across a lattice site. Taking D ∼ 10 μm²/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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Table of parameters and values used in cellular Potts model (CPM) simulations.

See text for more detailed reasoning behind values.

https://doi.org/10.1371/journal.pcbi.1006961.t001

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⁻¹ and η = 107 μm⁻¹ 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⁻¹. As seen in Fig 4E (gray), we find that the behavior of the DP depends sensitively on this timescale. When r⁻¹ is large, the DP increases with gradient strength. In contrast, when r⁻¹ 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⁻¹ = τ. 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 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⁻¹ 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⁻¹ is long, biased motion must be also persistent and vice versa. In contrast, when the memory timescale r⁻¹ 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 (8) where θ′ 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 θ′ and therefore promotes persistence. The normalization factor I₀ is the zeroth-order modified Bessel function of the first kind.

The requirement that P(θ|θ′) 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π)⁻¹, and non-negativity requires b < (2π)⁻¹ ≈ 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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Comparison of theory with experiments and simulations.

Colored circles show CI and DP for all values of bias parameter (color) and persistence parameter (size) for biased persistent random walk (BPRW) theory. Black lines show analytic approximations of the bounding curves. Red and cyan squares show experimental and simulation data, respectively, from Fig 4. Magenta numerals and arrows show “forbidden” regions and qualitative trends, respectively, discussed in text.

https://doi.org/10.1371/journal.pcbi.1006961.g006

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, (9) (10) respectively, where and z = χ/(1 − ψ), with and . We approximate the CI and DP in terms of these moments, (11) (12) 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 (left black line in Fig 6), where we have used the approximation M ≫ 1 (see Materials and methods). Note here that when CI = 0 we have DP ≈ (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 ∼ 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 in Eqs 11 and 12: because for any statistical quantity, we have , 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 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.

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 μ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 cm² T-flasks at 37 °C with 5% CO₂ 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% CO₂. 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⁻¹ 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⁻¹ 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: (13)

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 [56–58], so that the microfluidic platform could be maintained at 37 °C in a 5% CO₂ 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–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 (14) (15) where z = πb, and we have neglected the exponential terms in the limit M ≫ 1. Defining the small parameter ϵ = 1/M, these expressions become (16) (17) to first order in ϵ, where c ≡ 2(z⁻² − 2). Inserting these expressions into Eqs 11 and 12, we obtain (18) (19) to first order in ϵ. Because z and c are both functions only of b, we eliminate b from Eqs 18 and 19 to obtain (20) 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₀(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.