A CELLULAR POTTS MODEL SIMULATING CELL MIGRATION ON AND IN MATRIX ENVIRONMENTS

. Cell migration on and through extracellular matrix is fundamental in a wide variety of physiological and pathological phenomena, and is exploited in scaﬀold-based tissue engineering. Migration is regulated by a number of extracellular matrix- or cell-derived biophysical parameters, such as matrix ﬁber orientation, pore size, and elasticity, or cell deformation, proteolysis, and adhesion. We here present an extended Cellular Potts Model (CPM) able to qualitatively and quantitatively describe cell migration eﬃciencies and phenotypes both on two-dimensional substrates and within three-dimensional matri- ces, close to experimental evidence. As distinct features of our approach, cells are modeled as compartmentalized discrete objects, diﬀerentiated into nucleus and cytosolic region, while the extracellular matrix is composed of a ﬁbrous mesh and a homogeneous ﬂuid. Our model provides a strong correlation of the directionality of migration with the topological extracellular matrix distribution and a biphasic dependence of migration on the matrix structure, density, adhesion, and stiﬀness, and, moreover, simulates that cell locomotion in highly constrained ﬁbrillar obstacles requires the deformation of the cell’s nucleus and/or the activity of cell-derived proteolysis. In conclusion, we here propose a mathematical modeling approach that serves to characterize cell migration as a biological phenomenon in healthy and diseased tissues and in engineering applications.


1.
Introduction. Cell migration on and within tissues plays a critical role in a diverse array of processes, such as in developing embryos, where the coordinated movement of cells of different origin along extracellular matrix (ECM) layers is crucial for organogenesis, and migratory defects at all stages lead to severe embryonic malformations [42]. In adult organisms, cell movement is mostly quiescent, but active or re-activated, respectively, in processes such as immune surveillance or inflammation, where leukocytes migrate from blood vessels into infected tissues and then into the lymph node for effector functions [27], or in wound healing, where migration contributes to the repair of both basement membranes underlying epithelium and connective tissues. In pathologic conditions, cell migration is involved in chronic inflammatory diseases such as artheriosclerosis, or in cancer cell invasion into ECM and metastatization [69]. The process of cell migration is finally exploited in biomedical engineering applications for the regeneration of various tissues, such as cartilage, skin, or peripheral nerves in vivo or in vitro [11,38,75,86].
Cell motile behavior is determined by a complex set of mechanisms not only consisting of biochemical extracellular signaling and subsequent rearrangements, i.e. of the cytoskeleton, within the cell, but also of biophysical parameters of the surrounding environment, whose basic component is the Extracellular Matrix (ECM). The ECM forms a number of structures, i.e., two-dimensional (2D) flat basement membranes, or three-dimensional (3D) connective tissues that include 2D surfaces of big ECM bundles and can be described as a complex network of insoluble structural fibrous proteins such as collagen type I or others, glycoproteins and soluble glycosaminoglycans, which, together, provide microstructural guidance cues and biochemical stimuli for moving individuals.
How does a cell migrate then on and in ECMs? For the basic program of cell migration over flat ECM substrates, four requirements have been identified: migrating cells (i) are morphologically polarized in the direction of motion, (ii) adhere to their environment via adhesive molecules, i.e., integrins, (iii) generate traction forces by contraction of cytoskeletal elements, and (iv) retract their rear end [1,45]. For migration within 3D porous environments, in addition to these basic principles, the cell requires to steer its way through steric obstacles [17,66,83]. This can be achieved either by passing through constricted openings of the ECM by significant cell deformation and cytoskeletal force generation, or by activating a cell-derived proteolytic machinery able to degrade matrix components and to open space for cell movement [23,24,41,68,84].
This basic motile behavior is further modulated by a number of mechanisms that include determinants from both the surrounding extracellular matrix (adhesive ligands, fiber distribution, pore size, and elastic modulus) and the cell itself (adhesive strength, deformability, and proteolysis) (refer to [28,45,83] and references therein) that we aim to systematically analyze by a modeling approach based on an extended Cellular Potts Model (CPM, [4,33,34,36,50,71]). This is a grid-based Monte Carlo technique employing a stochastic energy minimization principle, used here to display the evolution of a cell population with distinct migratory behaviors that depend on matrix-or cell-derived parameters. As a distinct feature of our approach, each cell is modeled as a discrete compartmentalized object, differentiated into nucleus and cytosol, while the matrices are constituted of two components, an inhomogeneous fibrous collagen-like network, and a homogeneous interstitial medium.
The model is highly flexible, being able to characterize the migratory behavior of cells in several conditions, both on 2D substrata and in 3D ECMs. In the simulations, characteristics like cell shape and directionality are not imposed a priori, but are a result of the interaction with the matrix fibrous component. As an outcome, we focus on experimentally addressable characteristics of cell locomotion, i.e., cell overall displacement, velocity and persistence time, and cell shape, predicting how these quantities are influenced by manipulations of either cell or matrix properties.
Consistently with experimental observations, our findings provide evidence for a biphasic cell migratory behavior on planar substrates in response to variations of the number of matrix ligands or adhesion strength, with maximal movements at intermediate values. In 3D matrix environments, the geometrical distribution of the collagenous network, such as matrix alignment or pore size, or the matrix elasticity will be demonstrated to affect cell behavior in a similar way. Further, the cell compartmentalization allow to discern the mechanical rigidity of the nucleus that, being higher than the cytosol, limits the migration capacity of the entire cell. Finally, we will include ECM-directed proteolysis, resulting in enhanced migration in restricted environments.
The remaining parts of this publication are organized as follows: in Section 2 (Mathematical Model), we clarify the assumptions on which our approach is based. The computational findings are then presented in Section 3 (Results), where we separately analyze both topological and mechanical features of different matrices, and variations in cell biophysical properties on the cell migratory behavior. Finally, the results are discussed in Section 4 (Discussion), and in Section 5 (Appendix), we provide details on the estimates of statistical quantities used to characterize the migratory capacity of moving individuals.
2. Mathematical model. The above introduced cell-ECM system is modeled at the mesoscopic level using an extended Cellular Potts Model, a grid-based stochastic approach, which describes the behavior of single individuals and their interactions with the local microenvironment in energetic terms and constraints. The simulation domains are d-dimensional regular lattices (i.e., numerical repeated graphs formed by equivalent sites) Ω ⊂ R d , where d = 2, 3 (we will specify the spatial dimensions according to the specific application described in the following). Each d-dimensional site x ∈ Ω ⊂ R d is labeled by an integer number, σ(x), which can be interpreted as a degenerate spin originally coming from statistical physics [40,61]. As classically adopted in CPM applications, a neighbor of x is identified by x , and its overall neighborhood by Ω x , i.e. Ω x = {x ∈ Ω : x is a neighbor of x}. Subdomains of contiguous sites with identical spin form discrete objects Σ σ (i.e., Σ σ = {x ∈ Ω : σ(x) = σ}), which are characterized by an object type, τ (Σ σ ).
The spatial domain is then occupied by cells, ECM fibers and physiological liquid. The simulated cells, η, are defined as compartmentalized units, composed of two subregions which, in turn, are classical CPM objects Σ σ : the nucleus, a central cluster of type τ = N , and the surrounding cytosol of type τ = C. Each cell compartment is obviously characterized, as an additional attribute, by the cluster id η(Σ σ ) to identify the individual it belongs to. The cell population resides either on a 2D or within a 3D ECM.
The environment surrounding the cells is differentiated into a homogeneous medium-like state, τ = M , and an inhomogeneous collagen-like state, τ = F . The medium-like state reproduces the mixture of soluble components (among others, proteoglycans and glycoproteins in water), which compose the interstitial fluid constant in viscosity. The collagen state represents instead a network of insoluble macromolecules, that associates into first-order collagen fibrils and second-order fibers and displays the most abundant structure in mammalian tissues. Each fibrous component is treated as a CPM standard and non-compartmentalized CPM object Σ σ . Dimension, density and distribution of the fibrous structures will be specified in next sections and will reproduce 2D and 3D matrix types, respectively, typically employed for in vitro assays. The inclusion of an explicit two-component matrix environment, already present in some other CPM applications [5,32,49,67], is a fundamental aspect of this work: it allows an accurate analysis of how cell migratory behavior is influenced by the heterogeneous fibrillar extracellular environment and therefore by the ECM specific biophysical and biomechanical properties while cells glide in medium of constant and homogeneous physical properties.
The simulated cell culture evolves to iteratively and stochastically reduce the free energy of the overall system, defined by the so-called hamiltonian H, whose expression will be clarified below. The core algorithm is a modified Metropolis method for Monte Carlo-Boltzmann dynamics [36,52], which is able to implement the natural exploratory behavior of biological individuals. Procedurally, at each time step t in the model, called Monte Carlo Step (MCS), a lattice site, x source , is selected at random and assigns its spin, σ(x source ), to one of its unlike neighbors, x target ∈ Ω x : x target / ∈ Σ σ , also randomly selected. The net energy difference due to the proposed change of domain configuration, ∆H| σ(xsource)→σ(xtarget) = H (af ter spin copy) − H (bef ore spin copy) , is then evaluated. The trial lattice update is finally accepted with a Boltzmann-like probability function: where T Σ σ(xsource) (t) ∈ R + is a Boltzmann temperature. It does not reflect any conventional thermal temperature but it is a measure of the mobility of the moving compartment Σ σ(xsource) . The specific form of (1) is a definitive improvement of the classical function used in all CPM applications (formally recovered in the limit ε → ∞). The standard transition probability has in fact a significant weakness in the fact that, in the case of non positive net energy differences caused by the proposed displacement (∆H| σ(xsource)→σ(xtarget) ≤ 0), each element Σ σ(xsource) is certainly going to move, regardless of its intrinsic motility, given by T Σ σ(xsource) , which lacks biological realism. For example, a "frozen" cell (i.e., with negligible intrinsic motility) does not extend its pseudopods towards a chemical source even if it senses a high chemotactic gradient (which, in the absence of other external forces, would result in ∆H 0). This issue is addressed using transition probabilities similar to (1), i.e., which take into account the object motility T Σ σ(xsource) also in the case of energetically favorable displacement attempts. Indeed, the choice of function tanh is a modeling option: more in general, the reader can use any other continuous and increasing law p(T Σ σ(xsource) (t)) : as commented in detail in [71]. In particular, for τ (Σ σ(xsource) ) = N , T Σ σ(xsource) = T N gives a measure of the relative motility of the cell nucleus, while, for τ (Σ σ(xsource) ) = C, T Σ σ(xsource) = T C is a measure of the intrinsic motility of the overall individual, as it provides the frequency of the ruffles of its cytosol (which, on a molecular level, are determined by polarization/depolarization processes of the actin cytoskeleton, refer to [53,60,64] and references therein). Finally, for τ (Σ σ(xsource) ) = F , T Σ σ(xsource) = T F determines the vibration degree of matrix fibers. For each cell, T N is a low value (< 1), resulting in a more passive motion of the nucleus (with respect to the cell membrane), which, unable of autonomous movement, is dragged by the surrounding cytosol, characterized instead by a high T C 1 (see our recent work [70] for a detailed mechanical explanation). In most simulations, the matrix fibers are instead assumed to be fixed by setting T F = 0. For any given time t, the system hamiltonian, whose minimization drives the evolution of the system, is defined as: H shape models the geometrical attributes of simulated objects (both subcellular compartments and matrix threads), which are written as non-dimensional relative deformations in the following quadratic form: depending on the actual volume and surface of the object, v Σσ (t) and s Σσ (t) (which reduce, respectively, to its surface and perimeter in two dimensions), as well as on the same quantities in the relaxed state, V τ (Σσ) and S τ (Σσ) , corresponding to its initial measures. The formulation of (4) allows to have finite energetic contributions, as well as a blow up in the case of v Σσ (t), s Σσ (t) → 0, see again [71] for a detailed explanation. κ Σσ (t) and ν Σσ (t) ∈ R + are mechanical moduli in units of energy: in particular, κ Σσ (t) refers to volume changes, while ν Σσ (t) relates to the degree of deformability/elasticity of the related object, i.e., the ease with which it is able to remodel. Indeed, assuming that cells do not significantly grow or shrink during migration, the fluctuations of their volumes are kept negligible with high constant values κ Σσ = κ 1, for any individual η and for Σ σ such that τ (Σ σ ) = {N, C}. Moreover, cells moving in matrix environments are typically deformable, but their nuclei show a higher rigidity with respect to the cytoplasm region: therefore, for any η and for Σ σ such that τ (Σ σ ) = C, we set ν Σσ = ν C 1, while and for Σ σ such that τ (Σ σ ) = N , we set ν Σσ = ν N 1. The extracellular environment is instead assumed to have homogeneous mechanical and microstructural properties: in particular the matrix fibers are assumed to be rigid by setting κ F = ν F 1. However, it is useful to underline that in the following we will analyze how the explicit variation of fiber and nucleus stiffness will affect cell migratory phenotypes within 3D matrices.
H adhesion is the general extension of Steinberg's Differential Adhesion Hypothesis (DAH) [36,76,77]. In particular, it is differentiated into the contributions of either the generalized contact tension between the nucleus and the cytoplasm within the same cell, or the effective adhesion between a cell and both the medium and the fibrillar matrix component, and, in case of collision, between cells: The Js are binding energies per unit area, which are obviously symmetric. In particular, J int N,C implicitly models the forces exerted by intermediate and actin filaments, and microtubules to anchor the nucleus to the cell cytoskeleton and preventing cells from fragmenting. Even though cell-cell interaction is a very rare event, J ext C,C represents the local adhesive strength between neighboring cells, a measure of the local quantity of active and exposed cadherin molecules. J ext C,M and J ext C,F evaluate instead the heterophilic contact interactions between cells and matrix components.
Specifically, J ext C,M and J ext C,F are a measure of the affinity between cell surface adhesion complexes (i.e., sugar-binding receptors or integrins) to either non-solid (i.e., glycosaminoglycans in medium), or solid (i.e., fibrillar collagen) extracellular ligands, respectively [73]. In particular, given J int N,C 0 to prevent cell splitting, we assume J ext C,F < J ext C,M since, as widely demonstrated in literature, most cell lines in standard conditions adhere more strongly to the fibrous part of the extracellular matrix rather than to its soluble component (see [78] and references therein). J ext C,C is instead kept high to avoid cell-cell adhesive interactions upon accidental cell collisions that may affect cell movement. Setting constant and homogeneous values for the bond energies Js corresponds to the assumption of a uniform distribution of adhesion molecules on cell surfaces and of ligands in the external environment, without any change during the observation time. A summary of values of all the model parameters used in the simulations is given in Table 1.
Finally, it is useful to underline that, while in the 2D case cells can freely move on the entire extracellular ECM-coated surface, in 3D environments, the collagenous part of the matrix represents a potential steric obstacle that moving individuals must overcome during motion, whereas, in parallel, within interstitial medium they can freely float.

Simulation characteristics.
To apply the Cellular Potts Model to simulate and describe cell migration on and in ECM matrices, we start with default cell-ECM conditions and subsequently expose migrating cells to specific conditions, such as matrix orientation, density, adhesiveness etc., relevant for migration in vitro and in vivo.
The basic CPM in both 2D and 3D conditions contains certain common spatial and temporal characteristics. The spatial simulation domain Ω is a regular ddimensional lattice with periodic boundary conditions and a basic grid size (lateral edge) of 1.3 µm. In all the bidimensional simulations, Ω ⊂ R 2 represents a squared section, whose side length is 2.69 · 10 4 sites (i.e, 3.5 cm) of an experimental dish, which is commonly used for planar migration assays [16]. In the 3D case, Ω ⊂ R 3 reproduces instead an experimental scaffold with a volumetric extension of 1 cm 3 (formed by 4.55 · 10 11 cubic voxels). The temporal resolution of the model is the MCS which is set to correspond to 2 s to compare cellular dynamics with experimental observations. All the performed simulations last 12 h (≈ 21600 MCS) to ensure the development of sufficient long migration paths. The choice of this relatively short observation time allows also to exclude critical events, such as cell apoptosis and duplication, processes that ultimately go along with immobilization of cells [1]. As in particular during the approximately 1-3 hour-long mitotic phase cells undergo reversible rounding and abrogation of migration, inclusion of mitosis would reduce the sharpness of the readouts on migration dynamics and cell shape.
The basic cell-matrix model to be simulated consists of a heterogeneous ECM of fibrillar and amorph ('medium') components hosting a cell population of low density to allow isolated motions with very rare cell encounters.
In all 2D simulations, we plate 1 · 10 3 cells/cm 2 , as done in [16], while in all 3D simulations, we embed 2 · 10 3 cells/cm 3 , reproducing the cellular density of the experimental migration assays performed in [37]. The cells that interact with collagen-like fibers, i.e., fibroblasts or cancer cells of epithelial or mesenchymal origin, display initially a round non-migratory unpolarized morphology: therefore, as  default condition, we start with round flat disks with a central round nucleus for the modeling on 2D surfaces, and with spheres with a spherical nuclear compartment for the modeling in 3D matrices. In both cases, their overall diameter is 10 grid sites (≈ 14 µm), while the nucleus is 5 grid sites (≈ 7 µm) in diameter. For the reader's convenience, we underline that the entire volume (resp., area) of a cell is the sum of the volume (resp., area) of the nucleus and of the cytosolic region, while its external surface (resp., perimeter) is instead the difference between the surface (resp., perimeter) of the cytosol and the surface (resp., perimeter) of the nucleus. These dimensions, given in Table 1, reflect the mean measures of typical eukariotic cells except white blood cells [1]. In our model, we set the length of a collagen-like fiber equal to 15 lattice sites (≈ 20 µm). Its thickness would generally range between 100-500 nm [8,63,65], and therefore it would be substantially smaller than the grid resolution. However, following a common approach for CPM applications [32,67], we here accorded a fiber the measure of a single grid site, so that it is reproducible in the domain Ω. For sake of simplicity, we will use the term fiber for both the basic short ECM structure (≈ 20 µm long threads) simulated for the 2D condition, and the long structure crossing the entire spatial domain of the 3D cubic network. Depicted are standard two-component substrates containing both an isotropic fibrous ECM of moderate density (yellow stripes) and the medium (black), and cells (grey circles or spheres). As an initial condition, a sparse population of cells is plated on or into the matrices. Bottom panels (B): cell migration on or within the above-represented isotropic ECMs. Wind-rose graphs of 10 randomly chosen cell tracks over 12 h. Black circles represent the ending location of each cell center of mass. In both conditions, cells display a Brownian random movement with net final displacement ca. 50 µm, MSD ca. 9 · 10 4 ± 0.5 · 10 3 µm 2 (median 8.8 · 10 4 µm 2 ), and velocity ca. 10 ± 0.6 µm/h (median 9.7 µm/h). As reproduced from selected cell paths, the persistence time is low (ca. 1.5 ± 0.2 h, median 1.2 h). Here and in the following all values are given as means ± s.d. over 50 randomly chosen individuals (see appendix). The cell migratory behavior is consistent with the extracellular environment isotropy, and the absence of chemical gradients or other directional biases.

Isotropic 2D and 3D matrices.
We first test the model for standard matrices containing an isotropic, moderately dense, fibrous network for both two and three dimensions. As planar substrate, we distribute 3 · 10 5 flat collagen-like fibers in each x and y-direction of the 3.5 cm-side length dish, yielding a density of 500 fibers/mm 2 , see left top panel of Fig. 1(A). The analogous isotropic 3D scaffold consists of a regular cubic mesh of collagen fibers creating a uniform pore distribution of 10 µm width (i.e., the same order of magnitude of the initial cell diameter, see Fig.  1(A), right top panel). We simulate a regular fibrous network to avoid the minor heterogeneities often experienced in experimental matrices, where the distribution of the threads and the relative pore diameters is only roughly constant [47,54,82]. As shown in the wind-rose graph ( Fig. 1(B)), when cells migrate on both 2D and in 3D matrices, the selected cell paths display a random walk, without preferred direction, in the absence of biasing chemical gradients or matrix anisotropies.
Such migratory path structures and quantitative parameters are consistent with experimental results for both 2D and 3D porous ECMs for a large number of cell types from epithelial or mesenchymal origin. These include human adult vascular smooth muscle cells (HSMCs) plated on flat type IV collagen (CnIV) substrates of similar concentrations [16] or human glioma cells plated on polyacrylamide ECMs [79], and for different fibroblastic and cancerous cell lines migrating within 3D fibrous matrices of similar geometrical and structural properties, i.e., NR6 mouse fibroblasts in collagen-glycosaminoglycan matrices [37], or human melanoma cells in collagen lattices [25]. Indeed, these comparisons provide confidence in the choice of parameters describing the biophysical and mechanical properties of the simulated cell-ECM system.
3.3. Anisotropic 2D and 3D matrices. Next, we analyze the migratory characteristics of a cell population in the case of anisotropic matrices. In particular, we keep fixed the quantity of fibers as displayed in Fig. 1, but progressively change their distribution by increasing their number along the same x-direction, leaving remaining fibers disposed in their standard direction ( Fig. 2(A, B; top rows)). The alignment of the matrix is quantified by evaluating a proper index, that can be called alignment index, given by where d is the dimension of the domain, n x the number of threads along the xdirection and n tot their overall number. This quantity scales the percentage of fibers aligned along the x-direction, so that it is zero in the case of isotropic networks and 1 in the case of fully aligned matrices. As a result, for both 2D and 3D migration, the paths gradually adapt towards anisotropic random walks, in particular, the directional cell motion increases towards the principal direction of alignment, see Fig. 2(A,B; bottom rows). Interestingly, the cells final average velocity and MSD remain constant despite increasing substrate orientation, with very similar values for both 2D and 3D conditions (see Fig. 2(C)). However, the cell's 2D and 3D directed motile behavior in response to fiber distribution directly correlate with a strong increase in time (up to 5 hours) that cells are able to perform persistent (no back-and-forth) movement (see Fig.  2(C), right). Therefore, ECM geometry and architecture directly impact on the migration pattern of individual cells. The directionality of cell movement is here not introduced a priori, but is a direct result of well-defined directional-guidance cues provided by the specific matrices. The anisotropy of the matrices induces in fact a re-orientation of moving cells in the direction of the threads (i.e., with the formation of clearly distinguishable leading and trailing edges, see Fig. 2(A, B; top rows)) and the consequent motion along them, which is no longer an isotropic Brownian movement, but a highly biased locomotion.
The efficacy of cell migration is highly affected by the orientation and spacing of matrix components and its adhesive ligands, respectively, as experimentally proven by lithographic and microprinting techniques creating 1D ECM pathways that offer geometric guidance and adhesive structures at a microscale [9,17,19,46]. Moreover, several experimental models have demonstrated the cell's preference to migrate along aligned matrix fibers within 3D environments, such as fibroblasts in collagen [15] or neuronal cells in fibrin substrates [18]. Lastly, in vivo intravital imaging studies of carcinoma cells in the mammary fat pad have pointed out the preferential chemotactic movement of invasive malignant cells along thick bundles of collagen fibers offering a 2D surface towards blood vessels [13], while in the lymph node paracortex, the aligned microarchitecture of collagen and fibronectin fibers ensheathed by fibroblastic reticular cells significantly influenced the migratory behavior of T-cells [3].
3.4. Pore size in 3D matrices. The ECM fibers and bundles in in vivo tissues that provide guidance as well as barriers into structures that create pores and gaps of strongly varying local densities [82]. Connective tissues, i.e. of the skin, are categorized into loose and dense extracellular tissues and form irregular gaps between ca. 1 to 1000 µm exceeding, matching or measuring below a moving cell's diameter [82]. In diseased tissues, i.e., progressing tumors, the associated stroma often changes its architecture, i.e., into fibrotic tissues, over time [14]. Together, healthy as well as diseased tissues of varying gap sizes provide both guidance and physical barriers for moving cells.
In the CPM model, we simulate the effect of varying substrate fiber density on cell migration in 3D networks, where matrix fibers form a regular cubic mesh, with uniform pore sizes of originally 10 µm that increase of decrease due to modulation of the fiber numbers, whereas the diameter of the moving cell remains approximately 10 µm. In result, the simulations predict a bimodal behavior of cell velocity and persistence (Fig. 3(A)). At low numbers of fibers the 3D scaffold constitutes a sparse network, resulting in pores significantly larger than the diameter of the characteristic cell shape. In this case cells exhibit a short-range movement while their body remains in a stationary ameboid-like state, regardless of their deformation ability, presumably because the distance to the nearest distal matrix fiber is too high to experience adhesive interactions that enable cells to extend their membrane (Fig.  3(B)). On the other hand, the formation of pore diameters of cellular or slightly subcellular ranges allows cells to physically interact with fibers in all three spatial directions and is associated with most efficient migration rates [25]. In this case, migrating cells apply an elongated morphology and slightly reduce their diameters to ca. 8-10 µm. Finally, an increase in the abundance of 3D matrix threads results in the formation of a scaffold characterized by small pores with limited available space (i.e., half of a cell diameter or less), and a substantially decreased cell migration rate is predicted. The formation of long cytosolic formations are not sufficient to pass through such steric hindrance, as the nucleus can not enter (ν N it too high). A following section will examine the effect of nuclear deformability on cell migration. Cell elongation increases with decrements of pore size (i.e., increments of fiber number) until a threshold value, defined by the limit deformability of the nucleus. As in the following, each value in the plot is shown as mean ± s.d. over 50 randomly chosen individuals (see also appendix).
In summary, cells display a biphasic relationship that reveals most optimal migration at pore sizes at cellular or somewhat subcellular diameters, and diminishes at gaps greatly bigger or smaller than the moving cell diameter.
The outcomes of our models are consistent with the relative observations provided in the experimental literature. In 3D environments, neutrophil migration (both velocity and directional coefficient) has been reported to vary in a biphasic manner with the gel pore size [41], while mouse fibroblasts have been observed to migrate significantly more in collagen-glycosaminoglycan (CG) scaffolds featuring pore sizes slightly smaller than cellular dimensions, whereas they have exhibited less dispersion in matrices with larger pores [37].
3.5. Cell-Fiber adhesiveness. Cell-matrix adhesion is mainly mediated by integrins on the cell surface that connect ECM ligands to the cytoskeleton as well as to signaling molecules. Adhesion can be modulated by a number of parameters, such as (i) the number of substratum ligands, (ii) the expression and activation levels phology with fiber number on 2D substrates. The number of fibers is step-wise increased from 6 · 10 3 to 6 · 10 7 per dish (with 6 · 10 5 fibers per dish representing the standard case). All other parameters remain unchanged, such as in the standard case of Fig. 1. (A) The box-and-whisker plots (means are lines, medians are dots, see appendix) represent cell MSD, average velocity, and persistence time from 50 randomly selected cells. (B) Changes in aspect ratio (actual border length/circumpherence with same area) during migration over 12 hours upon varying fiber density. Number of fibers are: 6 · 10 3 (blue, low density), 2 · 10 6 (black, intermediate density), and 6 · 10 7 (green, high fiber density). Migration-associated lamellipodial ruffling is maximal at intermediate fiber densities, whereas at low and high number of threads cells remained roundish, associated with little migration. of integrins, and (iii) the resulting integrin-ligand binding affinity, which can be reduced by β1 integrin antibodies that block integrin binding epitopes to ECM or by soluble ligands that compete with ligand binding, or can be enhanced by integrin activating agents.
From the mathematical point of view, adhesiveness is modeled by both the substrate fiber density and the cell-fiber adhesion parameter J ext C,F .
3.5.1. Substrate Density of 2D Matrices. As mentioned above, adhesion depends on the number of substratum ligands applied here as varying fiber densities placed onto 2D surfaces. Since the variation of fiber concentrations in a 3D porous lattice will concomitantly change available space and thereby interdependent pro-migratory co-parameters, we excluded this approach from analysis. We simulate here both migration over a surface covered with an increasing amount of matrix fibers distributed equally and isotropically along the x− and y−directions, and related cell spreading. Cell spreading as an early adhesiondependent event is characterized by an increase of cell surface area over time [12], followed by the formation of migratory, i.e. pseudopodial structures. By keeping the variation of cell contact area with the underlying substrate nearly fixed (i.e., by high values of κ), our simulations do not capture cell spreading. Instead, we quantified elongation of pseudopodia as a measure for adhesion-dependent migration over 2D surfaces through the aspect ratio, defined in the appendix. Indeed, migration efficiencies develop a bell-shaped distribution from low towards high fiber numbers with a maximum at intermediate fiber numbers (Fig. 4(A)). At low ligand density, cells are unable to find sufficient collagen-like sites to attach and, in consequence, do not significantly displace. At the other extremum, an abundance of substratum ligands will lead to the formation of stable focal adhesions and, hence, low detachment and migration rates. Concomitantly, in both cases, cells mantain an ameboid-like shapes (see Fig. 4(B), lower inset). At intermediate fiber densities, relatively shortlived focal adhesions will form resulting in optimization of attachment-detachment cycles and maximal cell movement. The optimization of focal adhesion dynamics results in optimal cell movement and in an increment in membrane ruffling by the formation of membrane-rich structures, such as lamellipodia or filopodia, indicative of a migratory phenotype (see Fig. 4

(B), upper inset).
Different studies have coherently shown that migration on planar substrates with low fiber densities is limited by the cells impossibility to form sufficient attachments for the generation of traction and forward movement [24,35,44]. On the other hand, optimal ligand densities, precluding the formation of stable focal adhesions [2,12,43,57], cause rapid focal adhesion turnovers and result in maximal cell movement. Eventually, at high densities migration is blocked because integrin receptors engage into stable focal adhesions that exclude coordinated attachment-detachment for cell movement [21,30]. Blocked migration due to stable focal contact formation is usually accompanied by an increased spreaded area (again, refer to [12]), which we, however, did not capture with our approach.

Cell-Fiber Adhesion Strength for 2D and 3D
Matrices. Adhesion strength is mediated by integrin activation levels and resulting affinity to fibers. We thus simulate cell motility both over 2D surfaces and within 3D matrices (at standard conditions, Fig. 1) as a function of varying integrin-mediated cell-fiber adhesion strength J ext C,F . As a basic migration-adhesion relationship, the migratory capability of moving individuals can be sorted into the three regimes of high, intermediate or low adhesion strength and in principle is valid for movements both over a surface and within a 3D matrix (Fig. 5). At high integrin engagement (say J ext C,F < 3), cells display barely no detectable movement within the observation period, by being unable to detach from fibers. From an energetic view point, cells minimize the hamiltonian H by keeping such an adhesive contact. Given the high difference between J ext C,F and J ext C,M , moving individuals have in fact no benefit from further movements, meaning that an overly adhesive substrate causes the formation of integrin clusters on the cell surface strongly binding to substrate and not allowing detachment as needed for further migration.
Intermediate values of J ext C,F (say, in the interval [3,6]) yield moderately high adhesive forces associated with a balance of attachment and detachment, allowing cells to efficiently move along 2D surfaces or within the fibrous network with a maximal distance covered. Finally, above a certain value of J ext C,F (say, > 6.5), adhesion is lacking, and consequently cells (such as mesenchymal or epithelial cells) display barely detectable movement within the observation period. Given that J ext C,M < J ext C,F , cells prefer to fluctuate around the initial position in the interstitial fluid, avoiding contacts with the collagenous threads. Indeed, if a passive contact happens, cells soon detach from the fiber without exerting the traction needed for further movements.
The biphasic distribution of MSD is associated with a similar corresponding distribution in velocity, but only a flat curve in persistence (Fig. 5). Therefore, the adhesion-dependent overall motility is mediated mostly by a cell velocity, whereas the persistent component of cell motion consists of a relatively flat curve and refers with its low persistence level to random movement ( Fig. 1(B)).
In both two and three dimensions, the similarity of the biphasic dependence between the migratory properties of cells and their adhesiveness is consistent with published experimental literature, i.e., on tumor cells expressing high levels of β1 integrins [48]. However, they are not necessarily valid for all cell types, such as leukocytes that use adhesion-independent strategies when moving within a 3D collagen network [22]. The different assays used for cells when migrating on 2D or within 3D matrices have an impact on the conclusions of migration capacities. Whereas a non-adhesive cell detaches from a surface and cannot migrate anymore, non-adhesive cells are caught in a 3D network and may not, or may migrate by unspecific interactions with the lattice or by cytoskeleton-mediated propulsive mechanisms [24].
In accordance with our simulated data, a number of 3D ECM assays have shown similar trends for adherent cell types, such as human prostate carcinoma cells, whose velocity has been plotted as a biphasic function of an adhesiveness parameter such as ligands functionality as well as receptor density [88], or melanoma cells, cultured in collagen scaffolds and stimulated with different concentrations of integrin-binding peptide RGD [10]. Together, cell velocity can vary non-linearly with increasing ligand concentration, as it first increases, reaches a maximum and then decreases while the number of ligands still increase [47,57]. Thus, an intermediate level of cell adhesion to underlying or surrounding ECM is of crucial importance for the effectiveness of cell migration.  3.6. Fiber elasticity of 3D matrices. In the body, extracellular tissues display a range of elastic characteristics that are modulated by the collagen content, the amount of cross-links between collagenous molecules and the presence of elastic fibers. Dense tissues are usually rigid, thus increasing matrix density should enhance rigidity. In experimental studies using 3D ECM that where either modulated in density, i.e., fiber concentration [85], or rigidity [55,56,72,74], the other component becomes influenced as well. However, to separate at least theoretically these interdependent effects, we here simulate both varying scaffold stiffness (regulated by ν F ) and the geometrical microstructure. To quantify such convoluting factors, we provide contour plots of migration efficacy, as joint functions of pore size and fiber elasticity, that illustrate cell motile parameters as differently colored 'landscapes' (Fig. 6). It is useful to underline that elastic fibers are also characterized by a low constant T F = 0.2, as they are no longer rigid but can deform.
In Fig. 6(A) at high pore size (i.e., 20 µm), cells display a reduced motile behavior, regardless of the fiber stiffness, as already shown in Fig. 3(A). The rationale of this is that in very loose tissues cells migration is not supported by guiding fibers around the cell and consequently, the cell migrates along single fibers that, however, when stiff, promote to some extent traction and therefore migration. Next, at intermediate mesh dimensions, both cell velocity and persistence (and, consequently, overall displacement) reach maximal levels biphasically depending on matrix elasticity. If the collagenous threads are too elastic (i.e., ν F < 3), they can be easily deformed, without representing a sufficient anchor for pulling force generation required for cell motion. With a moderate stiffness (i.e., 5 < ν F < 9), the matrix fibers can be slightly arranged into contact guidance formations, thereby facilitating cell migration. On the contrary, a too rigid network (i.e., ν F > 9) forms steric obstacles that can be somewhat less efficiently overcome by moving individuals. Finally, small pore sizes allow motility only within elastic matrices, whereas migration is negligible for intermediate or high rigidities of the fibers. Migrating cells are able to move within small pores in fact only by significantly deforming the matrix network, creating open space to pass through. Therefore, if the pore size is much smaller than the dimension of the cell, matrix elasticity exerts an increasing influence. When evaluating the plots at constant rigidity, cell migration displays the same bimodal dependence on pore size previously captured in Fig. 3(A).
Such variations in fiber rigidity induce a suite of cell morphological changes (Fig.  7). Cells plated within rigid scaffolds typically display actomyosin traction-mediated elongated phenotypes, whereas, when cultured in progressively softer matrices they show instead decreasing elongation. Cells within complete compliant ECMs remain uniformly rounded, as cytoskeletal traction requires the counterforce provided by a rigid matrix [79]. In summary, cells migrate in a biphasic manner upon increasing separately density or stiffness, and also upon combined increase of density and stiffness reflecting experimental conditions best (imagine a decreasing curve in the plots of Fig. 6 from left top to right bottom).
As experimental examples, bimodal relationships between cell migratory ability and the deformability of 3D matrix scaffolds have been observed in experimental models of smooth muscle cells [59] and mouse fibroblasts, cultured in stepwise EDAC-cross-linked CG matrices of constant pore size [37]. A biphasic dependence on matrix rigidity has been previously reported also in isotropic homogeneous networks, as in the case of prostate cancer cells embedded in Matrigel with a fixed fibronectin level and variables stiffness [88]. 3.7. Nuclear compressibility in 3D migration. As pointed out in the previous section, to migrate within ECM of pores smaller than a cellular diameter, cells need to deform their body including the nucleus as the most rigid organelle [26]. Thus, the degree of nuclear deformability may contribute to the migration efficacy of a cell. The nucleus elasticity is mainly regulated by both the chromatin structure, and lamin intermediate filaments, that underlie the nuclear envelope [26,31]. The softness of a nucleus can be modeled by lowering the values of the nuclear rigidity ν N from 8.5 (see Table 1) to 0.5 (compare Fig. 6(A) and (B)). At high pore sizes of 10 µm (or higher) and low fiber rigidity (at left upper corner of left and middle panels in Figs. 6(A) and (B)), migration remains unaltered regardless of nuclear elasticity, as moving cells do not experience steric hindrance. However, once mesh dimensions and scaffold rigidity move to a intermediate ranges, nuclear elasticity somewhat facilitates cell movement, measured as MSD (an estimate product of persistence time and velocity). Whereas the velocity contributes to this enhanced overall movement, the persistence time remains equal (i.e., cells containing rigid nuclei already persistantly migrated through intermediate matrix networks). Lastly, as pore size and matrix elasticity further decrease to form a highly constrained environment (lower right corner), the simulations demonstrate that enhanced nuclear deformability provides a migration advantage characterised by both increased velocity and persistence time. Such a facilitated locomotion is mediated by an elongated and deformable nuclear configuration allowing the entire cell to squeeze and stretch more easily and thereby pass through the steric obstacles of a dense and rigid matrix (Fig. 8).
Our simulations relate to a number of experimental works, such as [81], where cell migration efficacy decreases with increasing matrix density and is associated with nuclear deformation, or [6], where glioma cell lines significantly deform their nucleus upon recruitment of non-muscle myosin II (NMMII) for squeezing through narrow locations in a brain model in vivo, thereby increasing their metastatic potential.
3.8. Matrix degradation in 3D migration. In the previous sections, we have demonstrated that cells move within matrix fibers of varying density and stiffness that act as constraints, deforming both body and nucleus. As an additional mechanism to overcome limited space cells may upregulate proteolytic enzymes that degrade ECM structure (i.e., matrix metalloproteinases, MMPs) that act either when bound to the cell surface or when secreted into the extracellular space. Accordingly, cells degrade steric fiber obstacles either in a cell contact-dependent manner targeting locally confining fibers, or, in a diffusive manner leading to gradient formation and consequently, a more overall weakening of the surrounding tissue structure [83]. As a result, barrier-free matrix spaces will be created resulting in longer distance traveling. We here perform simulations with cells that execute both contact-dependent and soluble proteolysis. The local concentration of the net protolytic activity (both surface-bound and diffusive) is defined as m(x, t), and is assumed to evolve following a standard reaction-diffusion equation: where δ(τ (Σ σ(x) ), M ) = 1 in the interstitial medium M and 0 elsewhere. λ m and D m are, respectively, the decay rate and the effective diffusion coefficient of proteolytic enzymes, constant and homogeneous in the extracellular environment. A low value of D m models proteolysis being strongly localized in regions close to cell membranes, in agreement with experimental evidence in [68,81]. P (x, t) models instead the local production of proteases either at the cell surface or secreted away from the external cell surface, at a constant rate π m which only activates upon local contact between the cell and a collagenous component: where we recall that C stands for cell cytosolic region. MMPs are capable to degrade the fibrous component of the matrix: to reproduce this biological effect, a lattice grid site x belonging to a degraded collagenous fiber becomes a generalized medium (fluid) site when its local level of MMPs (m(x, t)) is sufficiently high (in our simulations above 2.5 µM). This change is implemented by changing its type τ from F (fiber) to M (medium), as done in [32]. The comparison of cell migration of either MMP-active and MMP-inactive individuals (Fig. 6(A) and (C)), reveals that at high and intermediate pore size and/or low matrix rigidity, the proteolytic machinery does not appreciably affect cell motion. The loose fiber network does not represent a significative obstacle for cell migration, which therefore is not enhanced further by MMP activity. In the case of small pores formed by rigid collagenous fibers (lower right), MMP activity promotes instead cell migration appreciably. This suggests that proteases, by degrading matrix fibers, are able to break steric obstacles in the close proximity of moving individuals, opening spaces for them to sample greater distances without turning back. The role of MMPs activity in cell migratory behavior captured in our model is in good agreement with the experimental results provided in [62] for dermal fibroblasts embedded in molecularly engineered PEG hydrogels, where a significative increment in the number of migrating individuals was observed upon up-regulation of proteolytic enzymes.
In conclusion, summarizing all the examined parameters, cell migration is greatly influenced by a number of complex ECM-and cell-derived characteristics that, in addition, display a number of interdependencies [28] and, together, determine the net outcome on migration.
4. Discussion. Due to the increasingly recognized importance of cell migration processes in matrix environments and its exploitation for therapy and for tissue engineering, a growing number of theoretical models has been developed. These modeling approaches analyze the relative importance of single and interrelated parameters to predict migration behavior.
We employed a simple and intuitive version of the Cellular Potts Model to simulate the motile behavior of cells seeded either on 2D matrix substrates or embedded within 3D matrix scaffolds. In contrast to previous approaches, the Cellular Potts Model used here treats each cell as compartmentalized into nucleus and cytoplasm, whose movement is driven by explicit interactions with the extracellular environment differentiated into fibers and medium. The introduction of the nucleus and its mechanical properties on one side and of the extracellular matrix and its specific fibrous characteristics on the other side allowed to simulate for the first time both their specific contributions in cell migration.
In particular, we considered isolated pro-migratory parameters derived either from the ECM, such as orientation, pore size, ligand density, or rigidity, or from the cell, such as adhesion, nuclear rigidity, or proteolysis, that control both cell migration efficacy and migratory phenotypes. In all proposed cases, the computational results are consistent with a number of published experimental counterparts, and confirm linear or biphasic dependencies of migration dynamics on single matrix-or cell-derived determinants (see short paragraphs at the end of each result section).
As a clear advantage of a theoretical approach, we have been able to independently vary and modulate in a graded fashion all biophysical cell parameters and microstructural properties of the matrix environment, which is helpful in dissecting the complex relationships between cell motility and the biophysical, biochemical and molecular properties of the matrix [28]. However, a modeling approach that describes isolated parameters is unable to encompass the complexity apparent in migration processes in vivo. Some of these additional, here disregarded, factors are (i) additional matrix deposition of moving individuals, leading to altered traction generation, adhesion and contact guidance; (ii) soluble or matrix-bound gradients of chemoattractants; (iii) molecular signals transmitted from the ECM to cells (outside-in signaling), thereby changing the activity of polarization-or contractilitymediating proteins (Rac, Rho) [24]; or (iv) inside-out signaling for reinforcement of adhesion [79].
Despite the limitations of theoretical modeling, our approach could be applied to the design of synthetic implant materials, i.e., acellular scaffolds with optimal values of pore size and stiffness that may accelerate cell in-growth, critically for regenerative treatments [7,11,38,80]. Further, applying our modeling approach on defined cancer invasion and inhibition studies in vitro and in vivo may assist in predicting some outcomes on therapeutic interventions. At this regard, it would be biologically relevant to adapt our approach to specific cell lines, characterized by distinct biophysical phenotypes (i.e., intrinsic motility, elasticity, or proteases activity). This can be easily done by inheriting the model parameters from experimentallymeasured quantities, characteristic of the selected cell Further, it would be interesting to analyze collective migration of cellular aggregates fundamental in several physio-pathological processes, as commented in [39]. In such aggregates a differentiation may occur among individuals of the same origin, (i.e., tip and stalk cells during angiogenic processes, or leader and follower cells during a skin wound healing [29]), whereas competitions for nutrients or altered heterotypic interactions may significantly affect the migratory capacity of an entire cell lineage (for example, cancer cells of epithelial origin inhibit the motility and induce apoptosis in neighboring normal individuals). Obviously, in this case, it is necessary to define all cell types in the model framework, together with their phenotypic parameters and the mechanims underlying the behavior of and mutual interactions within the cell collective.
In summary, our findings may contribute to both understanding and exploitation of cell migration processes on and in tissues.
where t f inal corresponds to the final time of the observation period which, as explained in the text, is set to 21600 MCS (12 hours).
The mean squared displacement (MSD) of a cell η at time t is calculated as where x CM η (0) is the initial position of its center of mass. Following [16,87], the squared displacements are averaged over all previous time steps, in order to take into account the back and forth motions exhibited by the moving individuals. As demonstrated in a number of previous experimental [16,37] and computational [16,87] studies, at sufficiently long times the mean square displacements vary approximately linearly with the number of time steps. It can therefore be related to cell instantaneous velocity (v η ) and persistence time (p η , which quantifies the directional productive motion) with the so-called persistence-random-walk (PRW) law: In particular, at still longer observation periods, (12) reduces to: and the persistence time of a moving individual can be directly calculated as The PRW relation has been demonstrated to characterize the cells migratory behavior more properly than other common methods, which calculate the average distance migrated by biological individuals in an arbitrary time interval, as commented in [20]. For the statistical analysis, cells that do not display a final MSD greater than their diameters are classified as non-motile and assigned a velocity of 0 µm/h and an undefined persistence time, as we follow the criterion described in [16,30].
A.1. Statistics. Cell motile parameters (MSD, velocity and persistence time) are represented in the figures as box-and-whisker plots, where the edges of the boxes are the 25 th and 75 th percentiles and the whiskers the 10 th and 90 th percentiles. The horizontal line represents the median, while the large black dot corresponds to the mean of the distribution. Statistical significance (p < 0.05) was determined for motile fraction data by the Students' t-test and for non-normally distributed data sets by the Kolmogorov-Smirnov test over each 50 randomly chosen individuals. In the multidimensional contour plots the values of the cell migratory parameters are means over 50 randomly chosen individuals.
Quantitative analysis of cell morphological changes is carried out by evaluating the evolution of the cell aspect ratio, given by the ratio between the actual cell surface (respectively, border in 2D) and the surface of the sphere having the same volume (respectively, the length of the circumpherence having the same area in 2D). It is useful to underline that in our model cell volume (respectively, area in 2D) is kept nearly fixed by high values of κ in Eq. (4), see Table 1. Therefore the aspect ratio gives a quantitative measure of cell membrane ruffling. Finally, the time evolution of the aspect ratio is given in the plots with mean ± s.d. over 50 randomly chosen individuals.