Anomalous Diffusion as a Descriptive Model of Cell Migration

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 (MSD) 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 [1][2][3][4][5][6][7] and density of cell adhesive ligands [8][9][10][11][12] .
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 [13][14][15] , 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 2, [16][17][18][19] . 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. 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 20,21 .
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 [22][23][24][25][26][27] ; however, there are fewer examples of AD's use in describing adherent cell SANPAH (ProteoChem, Denver, CO; in pH 8.5 HEPES buffer) under UV light for 10 minutes, rinsed twice with PBS, and incubated overnight with 10µg cm -2 rat-tail collagen I.

Cell Migration
On coverslips, cells were seeded at 4000 cells per cm 2 and given 18 hours to adhere in growth medium. On 2D hydrogels, cells were seeded at 5,700 cells per cm 2 and given 24 hours to adhere in growth medium. In 3D hydrogels, cells were seeded at 1,000 cells µL −1 and allowed to equilibrate and spread within the gels for 24 hours, with a medium change to regular or conditioned medium 4 hours prior to microscopy. On coverslips with ECMcocktail 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 hours prior to microscopy. Brightfield and fluorescent images were taken at 15-minute intervals for 12 hours 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 minutes for 12 hours with 10 µm z-steps using a Zeiss Cell Observer Spinning Disc at 20X with a NA 0.5 air objective (Zeiss) and cells were tracked using Imaris. N ≥ 24 individual cell paths were generated for each condition, consisting of coordinates (x n (t), y n (t)) or (x n (t), y n (t), z n (t)) at observation times t = 0, Δt, 2Δt, … t max , where Δt is the time between image acquisitions, t max 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-hour observation period were excluded from all calculations. 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. Best-fit values of P (P n for cells and P agg for conditions) between 0 and t max were obtained using the Matlab (Mathworks) function lsqcurvefit(), which attempts to minimize the residual sum of squares (RSS), defined as

Displacement calculations of cell migration paths
The coefficient of determination (COD), a measure of goodness of fit, was calculated as is the total sum of squares (TSS) and the overbar (‾) indicates the unweighted mean for all τ (only up to 4 hours).
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 R 2 PRW . Generally, the use of the 2-parameter fit did not significantly increase R 2 PRW 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 (10) log $ AD ( ) = log Γ " Then log $ versus log data, up to a maximum τ of 4 hours, was fitted to the equation using lsqcurvefit() and 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 10000.
For AD model fitting of MSD vs. t 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 34 . The mean, median, quartiles, 10 th , and 90 th percentiles of individual-cell best-fit parameters were determined for each condition, and GraphPad Prism and Matlab were used to make the plots. Cellspecific parameters such as α n and P n were not weighted differently depending on the strength of the fit, i.e. R 2 value, when calculating condition-wide averages. Individual cells were considered subdiffusive for descriptive purposes if 0 ≤ α n < 1 and superdiffusive if 1 < α n ≤ 2.

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

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 2). Independently, we perturbed MDA-MB-231 chemokinesis and adhesivity, chemically, by adding either EGF (green) or a function-affecting antibody to b 1 integrin (red) (Fig. 2a-c).
On these rigid surfaces, regardless of the ECM protein cocktail or chemical perturbation, cells were largely (28-  (Table S1 and Fig. S1).
Across all conditions, both AD $ and PRW $ approached 1 as a n 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 AD $ and PRW $ remained greater than 0.95 for 97% of superdiffusive cells (a n >1), PRW $ decreased significantly as a n decreased below 1 (subdiffusive cells, Fig. S1). 82% of subdiffusive cells had AD $ >0.   Fig. S1).

The anomalous exponent, a, 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. 2). This is in line with previous studies of mammalian cell migration on other hard tractable ECM-coated surfaces 35 , which have previously been described with the PRW model. Cells exhibited much longer displacements on the lung ECM compared to brain and bone surfaces (Fig. 2a-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 b 1 -integrin targeting had a greater effect on the MSD, and on the AD model parameters α and Γ, than surface ECM.
The α n 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 ( Fig. 2d-f, 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-b 1 integrin to 84% on bone ECM-like surface with no chemical perturbation (average 63%;  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 Fig. 2d-f with S2). Average cell speed was more sensitive to chemical perturbation than substrate adhesivity (Fig. 2g). 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, α increases 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-b 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. from µm 2 /hr 1.6 to µm 2 /hr 2.0 ). When segregated in this way, Γ n , much like α n , was higher for EGFtreated cells than for untreated cells and anti-b 1 -treated cells; and higher for cells on lung-like than on brain-or bone-like ECM cocktail (Fig. S.3).

Γ, not a, is sensitive to substrate adhesivity
Given the known dependence of cell migration speed on the density and type of adhesive surface ligands 35, 36 , 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 RGD (Fig. 3a, c, e, and g) or full length ECM protein fibronectin (Fig. 3b, d, f, and h) using the same silane coverslip chemistry depicted in Figure 1 and applied in Fig. 2. To determine if the results from Figure 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 Fig. 3e-f), with a much narrower α n range than observed for 231s on the 2D tissue-specific ECM-mimicking surfaces.
For both adhesive ligands used, aggregate MSD increased with ligand concentration (Fig. 3a-b). Median cell speed ( Fig. 3c-d) and median Γ n (Fig. 3g-h and S4) increased with increasing fibronectin concentration, while median α n (Fig. 3e-f and S4 and Table S1) and median persistence time P n (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 a nor persistence time. Regardless, due to the largely superdiffusive cell populations for all conditions tested, average individual R 2 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 1, 37 , 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 PEG-PC hydrogel system, which has independent control over modulus and the density of adhesive ligands on the substrate surface 32 . We varied the substrate modulus from 1 to 64 kPa, and coupled type 1 collagen to the surface to make it cell adhesive (Fig. 4). The MSD of the breast cancer cells had a biphasic dependence on substrate stiffness, with maximum MSD occurring on 18 kPa gels (Fig. 4a), and similar to previous reports on biomaterial surfaces, we observed a biphasic dependence of cell migration speed as a function of gel stiffness (Fig. 4b) 1 . Somewhat in parallel with this biphasic response, α n and S n slightly positively correlated among cells on low-and high-modulus surfaces and slightly negatively correlated among cells on medium-modulus surfaces (Fig. S6).
Most cells migrating on these 2D soft gels were superdiffusive, with α values averaging between 1.3-1.4, and no observable relationship with gel modulus (Fig. 4c and S5b). Instead, there was a visible correlation between cell speed and Γ, and were both maximized on the 18 kPa surface (Fig. S2, Tables S2-3, and Fig. 4d). We conclude, as in Fig. 3, that the Γ parameter (acting as a transport diffusivity) is associated with cell speed, and the α parameter, (the trajectory anomality 38 is not. Furthermore, cell migration was superdiffusive regardless of gel stiffness (Fig. 4d,   S6a), or ligand density (Fig 3c,g), for the range of conditions we tested. For all samples, the average individual R 2 AD was at least 0.95 and the average individual R 2 PRW 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 [39][40][41] and measured cell motility in 3D as we exposed them to different pro-migratory chemical stimulations (conditioned medium from patient cell cultures). The gels were crosslinked with MMP-sensitive peptides, but the mesh size was orders of magnitude smaller than a cell (~20-25nm). 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. 5a), although speed across the population was not dramatically affected. Mean speeds of cell populations slightly increased as a function of stimulation (Table S.2), and individual cells were found to be faster in conditions with medium from either Patient 1 or Patient 3 compared to all cells in normal growth medium (Fig. 5b).
The Γ parameter had low values overall compared to the cells in 2D environments (Figs. [2][3][4][5] and was less affected by these medium conditions, revealing a stronger relationship with dimensionality compared to medium (Fig. 5d). As shown in Table S.1, the AD model fits exceeded the PRW fits for all conditions. AD and PRW fit equally well (around 0.95, Table S.1) 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 exception to this rule is best shown in the experiments from Figures 3 and 4, 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 were subdiffusive.
Comparisons between the 2D (Fig. 2) 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 often 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.
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 42 . For diffusive and superdiffusive cells, P varies positively with α, which was generally increased by EGF and by MSCconditioned 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 43 . 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. 2). This is not surprising given the chemo-kinetic potency of EGF 44 .
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). α increased significantly for cells in Patient average individual Γ increased for Patient 1 and Patient 3 MSC-conditioned medium compared to control medium (Table S3, Fig. S7). The AD parameters α n , Γ n , α agg , Γ agg similarly increased in EGF-treated medium and decreased in anti-β1 conditions (Table S1, S3; Fig. 2d-f, S3). Thus Γ and α are each independently responsive to chemical perturbation. On the other hand, α does not strongly depend on substrate adhesivity (Fig. 3e-f) or elasticity (Fig. 4c) while Γ does (Fig. 3h, 4d).
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. 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.
We saw that the average value of Γ, which typically followed closely with cell speed, increases with increasing ECM ligand density and is biphasic with respect to elasticity. If high ligand density and a mid-range/optimal elastic modulus help the cell attain maximum traction, perhaps Γ is a reflection of its magnitude.
It should be noted that we observed a very strong correlation between α and P for all conditions (Fig. S3).
Furthermore, Γ n strongly correlated with S n for cells with similar α n (data not shown). Together, these results suggest that the pairs (α n , Γ n ) and (P n , S n ) 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.
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.

Supplementary Material
See the supplementary material for complete data for cell motility parameters and model fitting data. Figure 1: 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 Γ. , and best-fit estimates for a n as histograms (d-f), 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 b 1 integrin (red).   given for each medium condition tested.