Anomalously diffusing and persistently migrating cells in 2D and 3D culture environments.

Appropriately chosen descriptive models of cell migration in biomaterials will allow researchers to characterize and ultimately predict the movement of cells in engineered systems for a variety of applications in tissue engineering. The persistent random walk (PRW) model accurately describes cell migration on two-dimensional (2D) substrates. However, this model inherently cannot describe subdiffusive cell movement, i.e., migration paths in which the root mean square displacement increases more slowly than the square root of the time interval. Subdiffusivity is a common characteristic of cells moving in confined environments, such as three-dimensional (3D) porous scaffolds, hydrogel networks, and in vivo tissues. We demonstrate that a generalized anomalous diffusion (AD) model, which uses a simple power law to relate the mean square displacement to time, more accurately captures individual cell migration paths across a range of engineered 2D and 3D environments than does the more commonly used PRW model. We used the AD model parameters to distinguish cell movement profiles on substrates with different chemokinetic factors, geometries (2D vs 3D), substrate adhesivities, and compliances. Although the two models performed with equal precision for superdiffusive cells, we suggest a simple AD model, in lieu of PRW, to describe cell trajectories in populations with a significant subdiffusive fraction, such as cells in confined, 3D environments.

INTRODUCTION

Cell migration is integral to a variety of physiological processes including organ development, tissue morphogenesis, wound healing, and immune response. A greater understanding of the motility effects of environmental cues can inform the design of biotechnologies such as movement-directing scaffolds. Research into the relationship between cell migration and cues from the cellular microenvironment increasingly takes advantage of the capability to manipulate properties such as the extracellular matrix (ECM) compliance and density of cell adhesive ligands.

Descriptive (i.e., empirical) models of migration dynamics facilitate analysis of microenvironment dependence in part by assigning parameters to characterize cells, individually and in aggregate. One of the most commonly used models for describing individual cell migration in 2D is the persistent random walk (PRW) model, whose mathematical formulation was originally developed as modified Brownian motion. Until recently, the migration of adherent cells has been explored almost exclusively on 2D surfaces, but is now investigated in 3D as well, partly due to the advent of bioengineered environments capable of encapsulating cells and more closely capturing in vivo conditions. Despite its success on 2D surfaces, cell migration is often not well described by the PRW model at any appreciably long time scale in confined 3D environments. Indeed, 9%–46% of low persistent (P < 50 min) motile cells in 3D collagen gels have been found to be poorly fitted by the PRW model (R2 < 0.6). Given the increasing use of 3D environments to study cell movement, there is a need for a model that can effectively describe individual cell movement in both 2D and 3D environments. Furthermore, there is a need to be able to predict migration model parameters that vary based on easily quantifiable and controllable extracellular conditions such as growth factor concentration, ECM ligand density and composition, and ECM compliance, all of which are known to have a significant effect on cell migration.

No simple descriptive model is commonly used to capture a broad range of types of cell motion, from highly constrained to ballistic. We therefore adapted the anomalous diffusion (AD) model for individual and aggregate cell migration. In contrast to normal (free) diffusion, in which the mean squared displacement grows linearly with the time interval τ, in anomalous diffusion, the mean squared displacement grows as a power, , of the time interval, where 0 < α < 2, by definition lending this model the flexibility to describe both sub- and superdiffusive motion. Variants of anomalous diffusion, in which α may be constant or τ-dependent, accurately describe a variety of physical and biological phenomena; however, there are fewer examples of AD's use in describing adherent cell migration in the literature, and it has not been used to systematically analyze individual cell trajectories to the best of our knowledge.

Given that many cells migrating in 3D are subdiffusive, we undertook to systematically characterize the trajectories of individual cells (and aggregate sample-wide migration) under various extracellular conditions using the AD model. We found that PRW and AD gave similar correlation coefficients for superdiffusive cells, but that the AD model was better at describing subdiffusive cells. The AD parameter α more clearly differentiated subdiffusive cells from each other than did the PRW parameter P (persistence time). The AD parameters as well as the PRW parameters were found to predictably vary with geometry, elastic modulus, ECM composition, and ECM ligand density. Therefore, we suggest the AD model is a more robust model of individual cell movement, particularly in constrained, 3D environments.

RESULTS

The AD model outperforms PRW in describing individual subdiffusive cell motion

We first quantified cell motility on supra-physiologically stiff surfaces: 2D coverslips coupled with full-length, integrin-binding (ECM) proteins. We created three different surfaces, inspired by proteins found in different tissues of the human body: bone, brain, and lung (Fig. 1). Independently, we perturbed MDA-MB-231 chemokinesis and adhesivity, chemically, by adding either epidermal growth factor (EGF; green) or a function-affecting antibody to β1 integrin (red) [Figs. 1(a)–1(c)]. On these rigid surfaces, regardless of the ECM protein cocktail or chemical perturbation, cells were largely (28%–84%) superdiffusive [1 < α < 2, Figs. 1(d)–1(f), Table S1]. Most cells were well described by the AD model, with individual  ≥ 0.8 for >90% of cells in this experiment and this study (Table S1 and Fig. S1).

FIG. 1.

The alpha parameter of the anomalous diffusion model is highly sensitive to chemokinetic perturbation. (a)–(c) The mean squared displacement is shown for all cells and time intervals on a log-log plot. For each graph, the solid line is the model fit for the anomalous diffusion model (AD), and the dashed line is the model fit for the persistent random walk model (PRW). (d)–(f) Best-fit estimates for αn as histograms, for the different ECM surfaces created: bone [(a), (d)], brain [(b), (e)], and lung [(c), (f)]. (g) Individual average cell speed box-and-whisker plots for all 9 substrate and treatment conditions. Control experiments (performed in standard growth medium) are shown in blue, supplemented with EGF (green), and treated with a function-affecting antibody to β1 integrin (red). N ≥ 88 cell paths were analyzed for each condition. Error bars in (a)–(c) are SEM. In g, boxes show 25th–75th percentile, and whiskers show 10th–90th percentile.

Across all conditions, both and approached 1 as α approached its maximum of 2 (Fig. S1). Given the flexibility of fitting for PRW, and that both models fit well, this is an argument for using PRW for cells on rigid 2D surfaces. While individual and remained greater than 0.95 for 97% of superdiffusive cells (α > 1), decreased significantly as α decreased below 1 (subdiffusive cells, Fig. S1). 82% of subdiffusive cells had > 0.8, while 45% of subdiffusive cells had > 0.8, indicating that the AD model fits cell mean square displacement (MSD) data far better than the PRW model. Therefore, when considering all cells from each population, the AD model outperformed the PRW model across all ECMs (Table S1). Finally, fitting with the AD model on an individual cell basis revealed population heterogeneity (in terms of the distribution of α and Γ within a population compared to the distribution of P and S) that persisted over the tracked time interval, which would have been missed by analysis with PRW alone (Tables S1–S3, Fig. S1).

The anomalous exponent, α, is most sensitive to integrin-binding to 2D rigid substrates

We observed sub- and super-diffusive behavior in individual cells during 2D migration on bone-, brain-, and lung-inspired-coated glass coverslips (Fig. 1). This is in line with previous studies of mammalian cell migration on other hard tractable ECM-coated surfaces, which have previously been described with the PRW model. Cells exhibited much longer displacements on the lung ECM compared to brain and bone surfaces [Figs. 1(a)–1(c)]. Blocking β1 integrin suppressed migration, while EGF treatment enhanced it on all surfaces, and lung-like surfaces tended to have greater migration-enhancing abilities than brain- and bone-like surfaces. When comparing the effects of these soluble factors versus the cell-adhesive proteins on the surfaces, EGF and β1-integrin targeting had a greater effect on the MSD, and on the AD model parameters α and Γ, than surface ECM.

The α distribution within each condition was typically unimodal and sensitive to the ECM adhesivity and soluble factors, highlighting the capability of the power-function model to describe a heterogeneous population of cells [Figs. 1(d)–1(f) and S2]. Regardless of the ECM protein cocktail or chemical perturbation, cells' individual anomalous exponents spanned the entire possible range 0–2 but tended to have a majority of superdiffusive cells, with superdiffusive fraction ranging from 28% on brain ECM-like surface with anti-β1 integrin to 84% on bone ECM-like surface with no chemical perturbation (average 63%; Table S1). Furthermore, all 9 conditions had an aggregate anomalous exponent greater than 1, indicating aggregate superdiffusive movement, with αagg ranging from 1.02 to 1.58 (average 1.40, Table S3).

The distribution of both αn and Γn shifted in response to soluble factor treatment, while αn was more sensitive to substrate variation than Γn [compare Figs. 1(d)–1(f) with Fig. S2]. Average cell speed was more sensitive to chemical perturbation than substrate adhesivity [Fig. 1(g)]. Treatment with anti-β1 integrin decreased median α, cell speed, and MSD across all three substrates, while EGF treatment increased median α and average displacement. Overall, for sample populations where both PRW and AD fit well, α increased with cell speed and persistence time, while there were no observed trends with Γ. Across all cell data, as expected, persistence time and α strongly positively correlated (Fig. S3), as did Γ and speed.

On average, a greater fraction of cells were superdiffusive in EGF-treated conditions (73%) than in untreated (69%) or anti-β1-treated conditions (48%). To fairly compare the effect on the anomalous diffusion coefficient Γn of different stimulants, we segregated individual cells by αn so that individual values of Γn would have approximately the same units (e.g., between μm2/h1.6 and μm2/h2.0). When segregated in this way, Γn, much like αn, was higher for EGF-treated cells than for untreated cells and anti-β1-treated cells, and higher for cells on lung-like than on brain- or bone-like ECM cocktail (Fig. S3).

We repeated this experiment with 231s and hTERTs on brain cocktail-coated glass coverslips for 24 to 48 h in order to understand to what extent the results hold for longer time interval ranges (τmax of 4 and 8 h). With both 4 h (16 MSD data points) and 8 h (32 MSD data points) fittings of the same cell trajectories, the average individual correlation coefficient for AD was higher than that for PRW in all the tested conditions (Table S6). Average R2 values for 8 h fittings were high enough for practical purposes but, as expected, were lower than those for 4 h fittings.

Γ, not α, is sensitive to substrate adhesivity

Given the known dependence of cell migration speed on the density and type of adhesive surface ligands, we examined the quality of the AD model, in comparison to the PRW model, on surfaces of varying densities of adhesive ligand. We created coverslip surfaces coupled with either the integrin-binding peptide Arg-Gly-Asp (RGD) [Figs. 2(a), 2(c), 2(e), and 2(g)] or full length ECM protein fibronectin [Figs. 2(b), 2(d), 2(f), and 2(h)] using the same silane coverslip chemistry depicted in Fig. 1 and applied in Fig. 2. To determine if the results from Fig. 2 were limited to the epithelial breast cancer cell line profiled, we expanded the study to a bone marrow-derived mesenchymal stem cell line (MSC). Virtually all cells were superdiffusive [Table S1 and Figs. 2(e)–2(f)], with a much narrower α range than observed for 231s on the 2D tissue-specific ECM-mimicking surfaces.

FIG. 2.

Gamma follows MSD and cell speeds during haptokinesis. Average MSD versus time intervals (log-log plots) for hTERT MSC cell migration on surfaces functionalized with different concentrations of RGD (a) or fibronectin (b). For each graph, the solid line is the model fit for the anomalous diffusion model (AD), and the dashed line is the model fit for the persistent random walk model (PRW). Individual cell speed box-and-whisker plots are shown for RGD (c) and fibronectin (d) surfaces. (e)–(h) The anomalous diffusion parameters αn (e)–(f) and Γn (g)–(h) are given for RGD (e), (g) and fibronectin (f), (h) surfaces. N ≥ 68 cell paths were analyzed for each condition. Error bars in (a) and (b) are SEM. In (c) and (d), boxes show 25th–75th percentile, and whiskers show 10th–90th percentile. Error bars in (e)–(h) are 95% confidence intervals.

For both adhesive ligands used, aggregate MSD increased with ligand concentration [Figs. 2(a) and 2(b)]. Median cell speed [Figs. 2(c) and 2(d)] and median Γn [Figs. 2(g)–2(h) and S4] increased with increasing fibronectin concentration, while median αn [Figs. 2(e) and 2(f) and S4 and Table S1] and median persistence time P (Table S2) showed no observable dependence on surface ligand concentration. Median αn was 1.73 to 1.76 for cells on RGD, and 1.71 to 1.73 for cells on fibronectin. These results suggest that greater ligand density increased random cell motility, and thus increased speed and Γ, with small changes to MSD. However, there was no effect on cell directionality, which explains why we observed no differences in α or persistence time. Regardless, due to the largely superdiffusive cell populations for all conditions tested, average individual R2 for both the PRW and AD models for all samples in this experiment was at least 0.98 (Table S1).

Γ and speed have a biphasic dependence on substrate stiffness

Given the known dependence of cell migration speed on substrate stiffness, we tested the AD model against cell migration data obtained from MDA-MB-231 breast cancer cells on substrates of varying stiffness. We used our previously published poly(ethylene glycol)-phosphorylcholine) (PEG-PC) hydrogel system, which has independent control over modulus and the density of adhesive ligands on the substrate surface. We varied the substrate modulus from 1 to 64 kPa, and coupled type 1 collagen to the surface to make it cell adhesive (Fig. 3). The MSD of the breast cancer cells had a biphasic dependence on substrate stiffness, with maximum MSD occurring on 18 kPa gels [Fig. 3(a)], and similar to previous reports on biomaterial surfaces, we observed a biphasic dependence of cell migration speed as a function of gel stiffness [Fig. 3(b)]. Somewhat in parallel with this biphasic response, α and S slightly positively correlated among cells on low- and high-modulus surfaces and slightly negatively correlated among cells on medium-modulus surfaces (Fig. S6).

FIG. 3.

Cell mean squared displacement, speed, and the diffusion coefficient are sensitive to substrate modulus. Average MSD versus time intervals (log-log plot) for breast cancer cell migration on hydrogel surfaces coupled with a constant density of collagen (10 μg cm−2) and varying substrate modulus (1–64 kPa). The solid line is the model fit for the anomalous diffusion model (AD), and the dashed line is the model fit for the persistent random walk model (PRW). Individual cell speed box-and-whisker plots are shown for each stiffness tested. (c) and (d) The anomalous diffusion parameters α (c) and Γ (d) are given for each stiffness tested. N ≥ 42 cell paths were analyzed for each condition. In (d), asterisks denote statistically significantly (p < 0.05) from value at 18 kPa. Error bars in (a) are SEM and in (d) are 95% confidence intervals. In (b) and (c), boxes show 25th–75th percentile, and whiskers show 10th–90th percentile.

Most cells migrating on these 2D soft gels were superdiffusive, with α values averaging between 1.3 and 1.4, and no observable relationship with gel modulus [Fig. 3(c) and S5(b)]. Instead, there was a visible correlation between cell speed and Γ, and were both maximized on the 18 kPa surface [Fig. S2, Tables S2 and S3, and Fig. 3(d)]. We conclude, as in Fig. 2, that the Γ parameter (acting as a transport diffusivity) is associated with cell speed, and the α parameter (the trajectory anomality) is not. Furthermore, cell migration was superdiffusive regardless of gel stiffness [Figs. 3(d) and S6(a)], or ligand density [Figs. 2(c) and 2(g)], for the range of conditions we tested. For all samples, the average individual R2AD was at least 0.95 and the average individual R2PRW was at least 0.77 (Table S1).

Cell migration in confined, 3D environments is largely subdiffusive

Cell motility models, such as PRW, were initially developed from data obtained from cells migrating on flat, 2D surfaces. However, in vivo, cell movement is largely in 3D, with cells surrounded by ECM and other cells. To test the effectiveness of the AD model in describing this type of confined cell movement, we used another PEG-based gel and measured cell motility in 3D as we exposed them to different promigratory chemical stimulations (conditioned medium from patient cell cultures). The gels were crosslinked with matrix metalloproteinase (MMP)-sensitive peptides, but the mesh size was orders of magnitude smaller than a cell (∼20–25 nm). Therefore, cell movement in this environment is confined, and cells are expected to exhibit subdiffusive motion.

We observed that all chemical stimulations increased the displacement of cells in 3D compared to normal growth medium [Fig. 5(a)] although speed across the population was not dramatically affected. Mean speeds of cell populations slightly increased as a function of stimulation (Table S2), and individual cells were found to be faster in conditions with medium from either patients 1 or 3 compared to all cells in normal growth medium [Fig. 5(b)]. 3D migration is largely subdiffusive (α < 1) in control medium and is approximately diffusive (α ≈ 1) in MSC-conditioned medium [Fig. 5(c)]. α increased significantly for cells in patients 1 and 2 MSC-conditioned medium compared to control medium.

FIG. 5.

Alpha describes subdiffusion of cells in soft, 3D hydrogels. (a) MSD as a function of time interval for breast cancer cells migrating in soft, 3D hydrogels supplemented with either normal growth medium (control), or conditioned medium (CM) from different sources of MSCs (Immortalized MSCs: IM, patient 1: P1, patient 2: P2, patient 3: P3). The solid line is the model fit for the anomalous diffusion model (AD), and the dashed line is the model fit for the persistent random walk model (PRW). (b) Individual cell speed box-and-whisker plots for the different medium conditions. (c) and (d) The anomalous diffusion parameters αn (c) and Γn (d) are given for each medium condition tested. N ≥ 24 cell paths were analyzed for each condition. Error bars in (a) are SEM. In (b)–(d), boxes show 25th–75th percentile, and whiskers show 10th–90th percentile.

The Γ parameter had low values overall compared to the cells in 2D environments (Figs. 1–3 and 5) and was less affected by these medium conditions, revealing a stronger relationship with dimensionality compared to medium [Figs. 5(d)]. As shown in Table S1, the AD model fits exceeded the PRW fits for all conditions. AD and PRW fit equally well (around 0.95, Table S1) for superdiffusive cells, but the goodness of fit for PRW rapidly dropped off as α decreased below 0.8 (data not shown). We conclude that the PRW model is not applicable for cells in 3D environments because of the very high proportion of subdiffusive cells. The AD model, as a power-function model, is the simplest adjustable monotonically increasing scale-invariant function, and is sufficient to capture the heterogeneity of these cell trajectories.

DISCUSSION

By analyzing different human cell lines in unique engineered biomaterial environments, we demonstrate that, overall, the AD model describes cell movement better than the commonly used PRW model. The superiority of the AD model compared to the PRW model was evaluated for each experimental condition using a student's t-test (Tables S4 and S6). Only minor differences were observed in the model fits on cells from the experiments from Figs. 2 and 3, which contained the highest proportion of superdiffusive cells. Superdiffusive motion is described equally well by both models, with PRW arguably giving more readily interpretable information (namely, persistence time and the relation between persistence time and adhesivity or elasticity). However, the AD model was equally good at relating motility behaviors to its own model parameters. The largest advantage of using the AD model was in 3D environments, wherein the largest population of cells was subdiffusive.

Comparisons between the 2D (Fig. 1) and 3D (Fig. 5) experiments highlighted the divergence in fitting between these two models because these were the conditions that contained the highest proportion of subdiffusive cells. PRW failed to describe the motion of cells in confined, 3D environments. In these cases, the AD model fit exceeded that of PRW. In 3D, where subdiffusive movement is expected, the PRW equation yields a poor fit, with P ≈ 0 such that either the trajectory is diffusive with speed S (the MSD function is linear) or the fitted parameters do not carry meaning. Cells in 3D confined environments have MSD functions unlimited to diffusive or superdiffusive paths, and thus need to be described by a model with flexible parameters, such as AD. Interestingly, the 3D gels also had the highest proportion of “diffusive” cells (α near a value of 1), as defined by 95% confidence intervals (supplementary material Table S5).

The effect of the environment on a cell's ability to move away from previous locations at a slower rate than if it “freely diffused” is reflected in the value of α. Although not explored here, high persistence and a strong correlation between cell speed and persistence are observed in both 2D free diffusion and for 1D-confined paths. For diffusive and superdiffusive cells, P varies positively with α, which was generally increased by EGF and by MSC-conditioned medium, and decreased by anti-β1 integrin antibody. Targeting integrin β1 increased the population of subdiffusive cells, supporting the known role for integrins as critical for polarization and persistent motion. The variation in the integrin-binding domains presented in different ECM cocktails also affected the MSD of cells, but treatment with EGF had the most significant impact on increasing the proportion of superdiffusive cells (Fig. 1). This is not surprising given the chemo-kinetic potency of EGF.

The range of α for cells in the 3D gels spanned from very subdiffusive (practically immobile) to very superdiffusive (ballistic motion, equivalent to moving in a straight line). Subdiffusive motion is typically associated with movement within a caged or confined environment. We speculate that tight mesh sizes of 3D hydrogels studied here recapitulate some aspects of a caged environment. Our 3D synthetic polymer-based gels have a relatively homogenous distribution of mesh sizes of approximately twenty nanometers (several orders of magnitude smaller than the smallest possible cell protrusion). Cells must degrade the gel to move in any direction. This reduces the speed and MSD as a function of time, and, at our imaging resolution, results in alpha values below that of free diffusion, or diffusion in a solvent. The AD model has been used in similar MMP-degradable PEG-based hydrogels in recent work, where cells exhibited anomalous diffusion with an exponent between 1.4 and 1.8. We observed cells doubling back along degraded tracks, leading to more subdiffusive behaviors than this recent report, but that α increased significantly for cells in patient 1 and patient 2 MSC-conditioned medium compared to control medium. Also, both aggregate and αn-grouped average individual Γ increased for patient 1 and patient 3 MSC-conditioned medium compared to control medium (Table S3, Fig. S7). The AD parameters αn, Γ, αagg, and Γagg similarly increased in EGF-treated medium and decreased in anti-β1 conditions [Tables S1 and S3; Figs. 1(d)–1(f) and S3]. Thus, Γ and α are each independently responsive to chemical perturbation. On the other hand, α does not strongly depend on substrate adhesivity [Figs. 2(e)–2(f)] or elasticity [Fig. 3(c)] while Γ does [Figs. 2(h) and 3(d)].

In this study, the AD MSD equation describes cell migration paths with excellent fit, and contains parameters that clearly depend on the biomaterial or growth factor condition. We compared the two models across a variety of cell culture environments, and with two distinct cell types (MSCs and MDA-MB-231s). It should be noted that these two cell lines will have different responses to adhesive ECM proteins, to growth factors, and to substrate stiffness, both in 2- and 3-dimensional environments. We argue that regardless of environment and cell type, compared to the PRW model parameters, the AD model parameters more accurately fit individual cell MSD data (Table S4) over 4 h and 8 h time interval fitting ranges (Table S6). While our study highlights the robustness of AD to describe cell movement, it is limited to using AD as a descriptive model and not a predictive one. In order to make a predictive model, future studies are needed to understand how α varies with experimental conditions, how Γ is tied to α, and how Γ varies with experimental conditions independently of α. The dependence of the AD parameters on substrate properties, growth factors, chemokines, and cytokines should be studied further to determine precise predictive correlations. Especially important are experiments with permeable 3D substrates to test the effects of varying elasticity and ligand density on subdiffusive movement. This would confirm that Γ is a strong proxy for cell motility in the absence of obstruction and not simply a reflection/artifact of high α.

A disadvantage of the AD model is that the parameters do not have readily apparent physical interpretations. Cell migration studies using the PRW model resulting in large differences in migration persistence time (and the measurable chemotaxis index) as a function of cell type or environmental perturbations could lead to mechanistic studies related to cell polarity and chemotaxis proteins. As a first pass of attempting to link the AD model parameters to environment, we observed that Γ often correlated with cell speed, and Γ increased with increasing ECM ligand density [Fig. 2(h)], and was biphasic with respect to elasticity [Fig. 3(d)]. This is a preliminary suggestion that Γ relates to cell traction forces or integrin-mediated signaling. Further studies could be performed to create better mechanistic links between Γ and α to cell-related processes and actionable hypotheses. A modification to the AD power law model via Lévy walk, although outside of this particular study, has been used in other cell tracking software, and could improve our understanding of how cells explore confined environments such as our 3D gel.

It should be noted that we observed a very strong correlation between α and P for all conditions (Fig. S3). Furthermore, Γ strongly correlated with S for cells with similar αn (data not shown). Together, these results suggest that the pairs (αn, Γ) and (Pn, S) contain similar information. Nevertheless, the AD model more conveniently represents data with a large fraction of subdiffusive cells because of the mathematical flexibility of α. The AD model thus has a broader “dynamic range” with respect to cell motility patterns and arguably reveals rather than masks heterogeneity within populations of subdiffusive cells. One additional notable difference between these two model approaches is that the PRW model is appropriate for all time scales, whereas the MSD vs tau log-log slope (i.e., alpha) can change with tau. This means piecewise AD fits are applied. If we look at how the MSD vs tau log-log slope changes with different time scales, we see a general trend that the alpha distributions shift to higher alphas for longer time observations (Fig. S9). This is an important consideration when designing experiments for tracking migrating cells and choosing time intervals between images.

Overall, our study highlighted the increased flexibility, and therefore better fit of the AD model compared to the more commonly used PRW model, particularly in instances for subdiffusive motion, such as in 3D environments. Overall we found that α was particularly sensitive to chemical stimulation (soluble factors in the medium or integrin inhibition), and Γ was more sensitive to substrate adhesivity and elasticity, tracking with cell speed for both. We recommend that the AD model can serve as a basis for simple computational models of cell speed in all environments, with specific applications in tumor growth, wound healing, or the colonization of an artificial tissue engineering scaffold seeded with cells.

METHODS

Cell culture

The human breast cancer cell line MDA-MB-231 used for cell migration on 2D coverslips and in 3D gels was a generous gift from Shannon Hughes, MIT. The MDA-MB-231 cells used for motility on 2D hydrogel surfaces were a generous gift from Sallie Schneider at the Pioneer Valley Life Sciences Institute. MDA-MB-231 cells were cultured in Dulbecco's modified Eagle medium supplemented with 10% fetal bovine serum (FBS), 1% penicillin-streptomycin (P/S), 1% l-glutamine, and 1% non-essential amino acids (NEAA, Thermo Fisher Scientific, Waltham, MA). Cells were cultured at 37 °C and 5% CO2. In the cases where MDA-MB-231 cells were cultured on 2D hydrogel substrates, l-glutamine and NEAA were omitted. Immortalized human mesenchymal stem cells (hTERT) were a generous gift from Linda Griffith, MIT, and were cultured in Dubelco's Modified Eagles Medium (DMEM) supplemented with 10% FBS, 1% P/S, 1% l-glutamine, 1% NEAA, and 1% sodium pyruvate. Patient-derived human mesenchymal stem cells were provided by the Texas A&M Institute for Regenerative Medicine and cultured between passage 1 and 5 in alpha Minimal Essential Medium (αMEM), 16.5% FBS, and 3 mM l-glutamine (Thermo). In cases where we created conditioned medium from other cell cultures, cells were seeded at 45 000 cells per cm2 in 2.5% αMEM with 3 mM l-glutamine for 72 h. The medium was removed and filtered through a 0.44 μm PES filter (Thermo) before use.

Preparation of ECM-coated coverslips

Coverslips were modified to present peptides and full-length ECM proteins as previously described [Fig. 4(a)]. Briefly, glass coverslips were silanized through vapor phase deposition of (3-aminopropyl)triethoxysilane (Sigma-Aldrich, St. Louis, MO, USA), rinsed, functionalized with 10 g L−1 N,N-disuccinimidyl carbonate (Sigma), rinsed, and functionalized through reactive amines with indicated amounts of RGD or fibronectin or ECM protein cocktails that were inspired by the ECM of tissues as follows: bone: 5 μg cm−2 of 99% collagen I and 1% osteopontin; brain: 1 μg cm−2 of 50% fibronectin, 25% vitronectin, 20% tenascin C, and 5% laminin; and lung: 2 μg cm−2 of 33% laminin, 33% collagen IV, 15% collagen I, 15% fibronectin, and 4% tenascin C (all weight%). Rat-tail collagen I and natural mouse laminin were purchased from Thermo; human tenascin C, human vitronectin, and human osteopontin from R&D Systems (Minneapolis, MN, USA); human collagen IV from Neuromics (Edina, MN, USA); and human plasma fibronectin from EMD Millipore (Billerica, MA, USA). Coverslips were then back-filled with PEG12 (Thermo) to block non-specific protein adsorption on any remaining surface area.

FIG. 4.

Tunable surfaces and modeling approaches to quantify and describe cell movement. (a) Overview of the chemistry used to create ECM-modified coverslips. A three-step process based on silane treatment results in protein- or peptide-modified surfaces (a generic peptide-modified surface is drawn). (b) Fluorescence results showing control of peptide surface coupling using the chemistry in (a). Results are from a model peptide (TAMRA-lysine) and read on a fluorometer. (c) Theoretical MSD plots for mildly (red), moderately (blue), and highly (green) persistent cell populations (S = speed, and P = persistence time) following the PRW model. (d) Theoretical MSD plots for subdiffusive (red), diffusive (blue), and superdiffusive (green) cells following the AD model, each with different, constant (time interval-independent) α and Γ.

A fluorescent assay [Fig. 4(b)] was used to quantify coupling of a generic, amine-containing peptide (TAMRA-lysine, Anaspec, Fremont, CA) to the coverslip surfaces, as a surrogate for peptides and the proteins used in the reaction depicted in Fig. 4. The concentration of TAMRA-lysine was titrated, and fluorescent intensity was measured in a fluorescent plate reader and measured at an excitation/emission of 545/575 nm.

2D PEG-PC hydrogels

Glass coverslips (no. 1.5 coverslip glass; Thermo) were UV-ozone treated (Bioforce Nanosciences, Salt Lake City, UT) and subsequently methacrylate-silanized with 2 vol. % 3-(trimethoxysilyl) propyl methacrylate (Sigma) in 95% ethanol (adjusted to pH 5.0 with glacial acetic acid) for 2 min, washed 3 times with 100% ethanol, and dried at 40 °C for 30 min. PEG-dimethacrylate (PEGDMA) (Mn 750, Sigma), from 0.6 to 9.1 wt. %, was combined with 17 wt. % 2-methacryloyloxyethyl phosphorylcholine (PC) (Sigma) in phosphate buffered saline (PBS). These PEGDMA crosslinker concentrations tune the Young's moduli of the resulting gels from 0.5 to 4.0 wt. % (1 to 64 kPa). Solutions were sterilized with a 0.2 μm syringe filter (Thermo) and degassed by nitrogen sparging for 30 s. Free-radical polymerization was induced by addition of 0.8 wt. % Irgacure (BASF, Florham Park, NJ). Hydrogels of 50 μl per coverslip were covered with an untreated coverslip and polymerized under UV light for 20 min. Post-polymerization, hydrogels were allowed to swell for 24 h in PBS, then treated with 0.3 mg ml−1 of sulfo-SANPAH (ProteoChem, Denver, CO; in pH 8.5 HEPES buffer) under UV light for 10 min, rinsed twice with PBS, and incubated overnight with 10 μg cm−2 rat-tail collagen I.

3D PEG-maleimide hydrogels

3D PEG-maleimide (PEG-Mal) hydrogels were prepared at 20 wt. % with a 20 K 4-arm PEG-Mal (Jenkem Technology, Plano, TX) and crosslinked in 2 mM triethanolamine (pH ∼ 7.4) at a 1:1 molar ratio with 50% 1.5 K linear PEG-dithiol (Jenkem), ∼17% GCRDVPMSMRGGDRCG, ∼17% GCRDSGESPAYYTADRCG, and ∼17% GCRDIPESLRAGDRCG (GenScript, Piscataway, NJ). Cells were pelleted and resuspended in PEG-Mal solution before mixing with PEG-dithiol. The cell adhesion peptides (GenScript, Piscataway, NJ) CGP(GPP)5GFOGER(GPP)5 (1 mM), CGPHSRN(G)6RGDS (0.8 mM), CGGSVVYGLR (0.1 mM), and CGGAEIDGIEL (0.1 mM) were reacted with PEG-Mal 10 min before gelation. Gels were polymerized in 10 μl volumes for 5 min before swelling in cell culture medium.

Cell migration

On coverslips, cells were seeded at 4000 cells per cm2 and given 18 h to adhere in growth medium. On 2D hydrogels, cells were seeded at 5700 cells per cm2 and given 24 h to adhere in growth medium. In 3D hydrogels, cells were seeded at 1000 cells μl−1 and allowed to equilibrate and spread within the gels for 24 h, with a medium change to regular or conditioned medium 4 h prior to microscopy. On coverslips with ECM-cocktail coatings, seeded MDA-MB-231 cells were treated with a live-cell fluorescent dye (CMFDA, Thermo) and then provided fresh medium or medium supplemented with 40 ng ml−1 epidermal growth factor (EGF, R&D Systems) or 0.83 μg ml−1 anti-β1 integrin (clone P5D2, R&D Systems) 4 h prior to microscopy. Brightfield and fluorescent images were taken at 15-min intervals for 12 h using an EC Plan-Neofluar 10 × 0.3 NA air objective (Carl Zeiss AG, Oberkochen, Germany) and cells were tracked using Imaris (Bitplane, St. Paul, MN, USA). Cell migration on 2D hydrogels was done in the absence of fluorescence, and cells were manually tracked. For 3D migration, fluorescent images were taken every 15 min for 12 h with 10 μm z-steps using a Zeiss Cell Observer Spinning Disc at 20× with a NA 0.5 air objective (Zeiss). Cell migration was quantified in using Imaris. Individual videos were drift corrected before analysis to mitigate drift, and background noise was measured by tracking a speck or bubble in the material to ensure only moving cells were considered for the analysis. N ≥ 24 individual cell paths were generated for each condition, consisting of coordinates [x(t), y(t)] or [x(t), y(t), z(t)] at observation times t = 0, Δt, 2Δt, … tmax, where Δt is the time between image acquisitions, tmax is the last acquisition time, and n ranges from 1 to N. Cells that contacted other cells, underwent division or apoptosis, or were not fully in frame for the entire 12- or 18-h observation period were excluded from all calculations.

Displacement calculations of cell migration paths

For reference, a description of all the pertinent symbols used in Eqs. (1)–(13) below is given in Table I. The individual mean square displacement (MSD) function of each cell was calculated for τ = 0.25 h, 0.50 h, 0.75 h, …, by averaging the square displacements over all available time intervals of length τ within the cell's trajectory (i.e., overlapping intervals) according to the equation where the angle brackets indicate averaging over all possible starting times t and the z-terms were neglected for 2D paths.

TABLE I.

Description of pertinent symbols used in the persistent random walk and anomolous diffusion calculations.

Symbol	Description
t	Time of observation
Δt	Time step between observations (e.g., 15 minutes)
tmax	Maximum value of t (e.g., 12 hours)
n	Cell index (e.g., 1, 2, 3, … )
N	Number of cells observed in a condition (maximum value of n)
(xn(t), yn(t) [, zn(t)])	Position coordinates of cell n at time t
τ	Time interval
d^2_n(τ)	MSD of cell n over time interval τ
d^2_agg(τ)	Aggregate MSD of cells in a condition over time interval τ
d^2^_PRWτ	Predicted MSD of the PRW model over time interval τ
P	Persistence time in the PRW model
Pn	Best-fit value of P for cell n
Pagg	Best-fit value of P for all cells in a condition
S	Root mean square speed
Sn	Root mean square speed of cell n (calculated or fitted)
Sagg	Root mean square speed of all cells in a condition (calculated or fitted)
RSS_PRW	Residual sum of squares for the PRW fit of the MSD as a function of τ
d^2(τ)¯	Mean MSD over all included time intervals τ (e.g., from 0.25 to 4 hours)
TSS_PRW	Total sum of squares for the PRW fit of the MSD as a function of τ
RPRW2	Coefficient of determination of the PRW fit of the MSD as a function of τ
d^2^_AD(τ)	Predicted MSD of the AD model over time interval τ
Γ	Non-linear diffusivity coefficient in the AD model
α	Anomalous exponent in the AD model
αn	Best-fit value of α for cell n (anomality of the pathway)
αagg	Best-fit value of α for all cells in a condition in aggregate
Γn	Best-fit value of Γ for cell n
Γagg	Best-fit value of Γ for all cells in a condition in aggregate
RSS_AD	Residual sum of squares for the AD fit of the MSD as a function of τ
logd^2(τ)¯	Mean of the logarithms of MSD over all included time intervals τ
TSS_AD	Total sum of squares for the AD fit of the MSD as a function of τ
RAD2	Coefficient of determination of the AD fit of the MSD as a function of τ

The aggregate MSD function of each condition was calculated according to where the angle brackets indicate averaging all available intervals of length τ from all available cell trajectories in the condition and the z-terms were neglected for 2D paths.

Fitting of or to the AD or PRW model functions was performed on points with τ values between Δt = 0.25 h and τmax = 4 h—only up to 1/3 of the maximum, so that each value of used in the fitting was obtained by averaging data containing at least 3 statistically independent displacement values. For model fittings with higher τmax, we also observed cells for between 24 and 48 h and performed AD and PRW model fittings with τmax of 4 h (16 lowest data points) and 8 h (32 lowest data points) on the resulting τ–MSD data. This experiment was done with the pooled data from 2 biological replicates, each with 2 technical replicates.

The persistent random walk model

Individual and aggregate MSD data up to a maximum τ of 4 h were fitted to the persistent random walk (PRW) equation used by Dunn in which S is root mean square speed over all available Δt-sized intervals (15 min) and P is persistence time, a characteristic parameter describing the average length of time over which cells resist major changes in direction of movement. For individual cells, the persistence is P and the root mean square speed is calculated as while for cell samples within a given condition, the PRW parameters are Pagg and In Eq. (3), S is the instantaneous speed, which could be different than the average speed calculated at t = 15 min in Eq. (4). We tracked cells and compared the root mean squared speed to the calculated speed at three time intervals (1, 15, and 60 min, Fig. S8). Longer time intervals resulted in closer correlations, and all the data in the main manuscript are at a time interval of 15 min. Best-fit values of P (P for cells and Pagg for conditions) between 0 and tmax were obtained using the Matlab (Mathworks) function lsqcurvefit(), which attempts to minimize the residual sum of squares (RSS), defined as The coefficient of determination (COD), a measure of goodness of fit, was calculated as where is the total sum of squares (TSS) and the overbar (−) indicates the unweighted mean for all τ (only up to 4 h).

Whereas the PRW model has only 1 fitted parameter (P), the AD model has 2. Since additional “degrees of freedom” can be expected to increase the correlation coefficient, some comparisons of the accuracy of the AD and PRW model might be reversed in favor of PRW if both S and P are allowed to vary. So a 2-parameter PRW fit was performed to obtain additional, arguably more “fair” values of S, P, and R2PRW. Generally, the use of the 2-parameter fit did not significantly increase R2PRW but sometimes drastically altered the best-fit parameters (see Table S2). Unless otherwise indicated, reported P values are from the 1-parameter fit.

The anomalous diffusion model

The anomalous diffusion (AD) model equation where is the anomalous exponent and Γ is the “anomalous diffusivity” parameter, was linearized by taking the logarithm to obtain Then data, up to a maximum τ of 4 hours, was fitted to the equation using lsqcurvefit() and Microsoft Excel's slope() and intercept() functions. AD parameters αn and Γ, for individual-cell MSD data, or α and Γagg, for aggregate MSD data, were simultaneously determined. Best-fit α was restricted to between 0 and 2 and best-fit Γ restricted to between 0 and 10 000.

For AD model fitting of MSD vs. τ data (both individual and aggregate), RSS, TSS, and COD were defined, respectively, as

Data and model analysis

A suite of purpose-built Matlab programs was used to process and analyse the migration data at the cell and aggregate level, including model fittings, and report, classify, and compare the results, and produce analytical plots such as Γ-binned-by-α plots and histograms. The MSD calculations and PRW model fittings were confirmed using the Visual Basic program DiPer. The mean, median, quartiles, 10th, and 90th percentiles of individual-cell best-fit parameters were determined for each condition, and GraphPad Prism and Matlab were used to make the plots. Cell-specific parameters such as α and P were not weighted differently depending on the strength of the fit, i.e., R2 value, when calculating condition-wide averages. Individual cells were considered subdiffusive for descriptive purposes if 0 ≤ α < 1 and superdiffusive if 1 < α ≤ 2.

Ethics approval

Ethics approval is not required for this study. Conditioned medium from patient cells described in Fig. 5 were de-identified and exempt from IRB approval.

See supplementary material for complete data for cell motility parameters and model fitting data.