Cell migration is a complex process involving many intracellular and extracellular factors, with different cell types adopting sometimes strikingly different morphologies. Modeling realistically behaving cells in tissues is computationally challenging because it implies dealing with multiple levels of complexity. We extend the Cellular Potts Model with an actin-inspired feedback mechanism that allows small stochastic cell rufflings to expand to cell protrusions. This simple phenomenological model produces realistically crawling and deforming amoeboid cells, and gliding half-moon shaped keratocyte-like cells. Both cell types can migrate randomly or follow directional cues. They can squeeze in between other cells in densely populated environments or migrate collectively. The model is computationally light, which allows the study of large, dense and heterogeneous tissues containing cells with realistic shapes and migratory properties.
Single migrating cells fall into different morphological categories [1]. For example, keratocyte-like cells, named after fish epithelial cells, retain a half-moon shape and glide quite persistently, led by a lasting, thin and broad lammellipodium [2,3]. In contrast, amoeboid cells change their shape frequently and migrate erratically by extending and retracting protrusions [4]. Leukocytes are classical examples of amoeboid cells patrolling through dense tissue. For instance, T lymphocytes squeeze through the tightly packed epidermis to find and fight infections [5–7]. Theoretical models are useful tools for exploring the relation between cell shape, cell migration and their tissue level implications, because they allow experiments on virtual cells in controlled environments. Categorized by scale, cell-based models can be single-cell or multicellular. Single-cell models answer questions about cell behavior by implementing various intracellular mechanisms, while multicellular models answer questions at the tissue level by using simplified representations of individual cells in order to reduce computational complexity. Several single-cell models explain the shape and behavior of individual cells biomechanically [2,3, 8–12] or by other intracellular mechanisms [13,14]. Although they vary greatly in their goal and level of detail, these models generally obtain migration by using a variant of the “local activation global inhibition” mechanism [15–18], and they usually require many equations and parameters in order to impose physical forces [2,3], to allow forces to emerge [8], and/or to diffuse intracellular molecules [13,14]. Multicellular models typically use simplified cell shapes—points, spheroids, collections of lattice sites—and impose cell migration by means of vectors or other directional signals. The few multicellular models that combine cell shapes with the associated migration patterns demonstrated that realistic cell representations can be essential for the behavior that emerges at the tissue level [7,19,20]. In general, single-cell models that result in realistic shape and migration are computationally expensive and therefore difficult to use in a multicellular context. Multicellular models, on the other hand, tend to strip the cells of their inherent shape-migration interconnection. We propose a simple and computationally light phenomenological model of cell migration, implemented within the Cellular Potts Model (CPM), that faithfully reproduces cell shapes and their associated migration properties, and can readily be used in big multicellular simulations. The cells in our model reproduce amoeboid and keratocyte behavior, migrate randomly with different levels of persistence, and display qualitative traits of chemotaxing cells. We illustrate the capabilities of our model in a multicellular context with 2 experiments: one experiment with T lymphocytes migrating in tightly packed skin, and one experiment with keratocytes migrating collectively. We show that our mechanism merely increases the computation time of the CPM by a small percentage, and that this percentage stays constant with increasing cell numbers, making the mechanism scalable and therefore suited for large simulations.
Our new computational model, the Act model, is inspired by actin dynamics. It extends the CPM with a local feedback mechanism resulting in cell protrusions and, as a consequence, in cell motility. The mechanism amplifies the inherent membrane fluctuations of CPM cells in a manner depending on the size of the fluctuations and their recent protrusive activity. The protrusive activity is tracked by keeping an activity value for every lattice site. The empty lattice sites that form the medium have a zero activity value, while sites that are freshly incorporated by a cell get the maximum activity value (MaxAct). The activity value of a site decreases by one after every MCS, until it reaches zero, creating a memory of MaxAct MCSs in which the site “remembers” that it was active. The combination of the memory of a site u with the activity in its neighborhood forms the basis for a local positive feedback mechanism that biases the copy attempt from the active site u to a less active site v. Formally, the mechanism is implemented by subtracting the term Δ H Act (u → v) = λ Act / Max Act (GM Act (u) - GM Act (v) ) from the energy difference of the system, Δ H. The parameter λAct is the maximum contribution of the Act model to the Hamiltonian. The average activity around site u is calculated as the geometric mean of the activity values in the neighborhood of u, GM Act (u) = (∏ y ∈ V (u) Act (y) ) 1 / | V (u) |, with V (u) the direct Moore neighborhood of u that belongs to the same cell as u and 0 ≤ GMAct (u) ≤ MaxAct. Compared to the arithmetic mean, the geometric mean favors the contribution to the Hamiltonian of complete Moore neighborhoods with consistently high activity values; it diminishes the contribution of Moore neighborhoods with low activity values and it nullifies the contribution of Moore neighborhoods “with holes” (i. e. , lattice sites with activity value zero). The dynamics of a copy trial are shown in Fig 1. ΔℋAct can be interpreted as the force resulting from pushing and resistance at the membrane element between u and v. This force reaches a maximum when u, backed up by a recent active neighborhood (GMAct (u) = MaxAct) tries to extend into a lattice site that opposes no resistance derived from activity (GMAct (v) = 0). The Act model encourages relatively big patches of active lattice sites just behind the cell membrane to extend further, creating protrusions in a way that roughly resembles protrusion formation in real cells, where branched networks of actin filaments grow behind the membrane and push it outwards. Within an actin network, the actin globular subunits attach at the barbed ends of the filaments, age and detach some distance behind the leading edge [21]. In our model, newly incorporated lattice sites also “age” and “detach” as their activity values decay until they reach zero and the sites cease to be part of the active patch. From this perspective, the MaxAct parameter can be interpreted as the expected lifetime of an actin globular subunit within an actin network and the λAct parameter can be seen as the maximum protrusive force of the network. The parallelism between actin globular subunits and active lattice sites motivates also our choice for the geometric average instead of the arithmetic average in the calculation of GMAct (u). New actin subunits contribute to actin network growth by attaching to the existing network structure, therefore, in our model, new active sites need to adjoin the active patch they originate from. Because the arithmetic mean disperses the activity values, it potentially creates active sites that are disconnected from active patches, while the geometric mean prevents this effect. Adding a new mechanism to an existing model comes at a computational cost. To assess both theoretically and empirically the effect of the Act model on the performance of the original CPM, we evaluated the asymptotic computational complexity of the system, and measured CPU time for several simulations with and without the Act model. One strength of the CPM is its capability to handle all cell properties locally, at the level of single copy attempts from one lattice site to another. In a CPM simulation using a 2D square lattice of width ℓ, ℓ2 copy attempts of lattice size independent computation time are performed during each MCS, hence, the asymptotic computational complexity of one MCS is O (ℓ 2). Because the Act model is also entirely based on local information, adding it to the CPM preserves the asymptotic complexity of the system, increasing only the time required to evaluate one copy attempt. In practice, the computation time of one MCS changes not only with the size of the lattice, but also with number of cells on the lattice. CPM implementations are optimized by discarding or completely avoiding trivial copy attempts within the same cell or within the medium, practically acting only at cell borders. Adding cells on the lattice increases the number of cell borders where the non-trivial and computationally more expensive copy attempts take place. Therefore, increasing the number of cells leads to an increase in the computational cost of one MCS. We empirically evaluated the cost of the Act model in such a setting by measuring the CPU time of the CPM with and without the Act model for an increasing number of cells on a fixed size lattice (Fig 12). In these simulations, the Act model increased the CPU time of the system by about 10%, and this overhead did not change appreciably with the number of cells. In conclusion, enriching the CPM with the Act model preserves its asymptotic computational complexity and increases its computational cost merely by a small percentage. Importantly, this percentage remains constant for large multicellular simulations. Within the widely used CPM framework, we have created the Act model based on an actin-inspired mechanism. In the new model, cell shape and motility emerge from the interplay between stochasticity, local positive feedback, and membrane tension. Act model cells break symmetry randomly and move in an amoeboid or keratocyte-like manner. When placed in a chemokine gradient, the in silico amoeboid cells reproduce several qualitative features of real chemotaxing amoeboids: they form more protrusions up the chemokine gradient and they react proportionally to the steepness of the gradient. In a multicellular context, in silico keratocyte-like cells reproduce the density dependent collective migration of real keratocytes [35], while amoeboid cells are able to squeeze themselves through densely packed tissues, resembling effector T cells crawling through the skin epidermis [7]. The Act model is computationally light, with its computational complexity scaling well with the number of cells on the lattice. Our model is conceptually related to existing single-cell models that reproduce amoeboid and keratocyte-like cell types [8,14]. Within a phenomenological framework similar to the CPM, Nishimura et al. [14] obtained cell migration by using the interaction between an actin surrogate and a generic actin inhibitor. In their model, cell deformations create regions of low concentrations of inhibitor which trigger actin surrogate accumulation and further deformations. This results in a feedback mechanism that drives cell polarization and movement. The Act model reproduces the same rich qualitative behavior with a much simpler implementation of actin-like dynamics—it uses two extra CPM parameters versus the six migration related parameters in the Nishimura model—and is sufficiently light to be used in various multicellular contexts. Multicellular models typically simplify cell shape and cell migration mechanisms in order to reduce computational complexity. For instance, in a 3D CPM model of T cell migration in lymph nodes, Beltman et al. [36] implemented cell movement by steering cells along direction vectors that are periodically updated to the recent cell displacement. This mechanism also results in persistent random migration, and has a similar computational complexity to that of the Act model, however, it cannot reproduce amoeboid behavior because the cells remain roundish. Importantly, in their model cell migration is imposed, while in the Act model cell migration emerges from the changes in cell shape. A limitation of the Act model is its phenomenological character, which allows only a simplified representation of the actin network and its dynamics. The Act model is implemented within the CPM, which is inherently phenomenological, although there are ongoing efforts to relate the parameters of CPM to biophysical properties of cells [32]. The forces acting upon CPM cells are implicit, controlled globally by the Hamiltonian, and cannot be imposed explicitly at the level of the cell. This means that most cell features can only be inferred by measurements performed during/after simulations. Therefore, the calibration of in silico cells amounts to finding the set of parameter values that results in the desired shape, biophysical, and/or migratory cell properties, with prior knowledge on how each CPM parameter influences certain cell properties (see also [37], part II). Noteworthy, the migration features of in silico cells in a non-empty environment are partly determined by the environment and therefore they typically cannot be set a priori, but need to be evaluated from simulations. This holds true for all phenomenological migration models proposed thus far. In the present implementation, the Act model resides on a 2D grid, making it especially suitable for simulating migration on flat surfaces or in very thin tissues, like epithelia. Adapting the model to 3D is straightforward and would facilitate the study of different 3D phenomena like amoeboid migration in soil, the individual or collective migration of cancer cells during metastasis, and lymphocyte migration in organs and tumors. In conclusion, our simple model is actin-inspired, achieves emergent amoeboid and keratocyte-like behavior with only two additional parameters and reproduces directed and persistent random migration. Because the Act model is light and simple, it can readily be used in large and heterogeneous multicellular simulations, while preserving the shape and behavior of the single-cell.
The Cellular Potts Model (CPM) [38] is a powerful lattice-based method for modeling cells and tissue dynamics, while retaining individual cell identity. CPM cells occupy areas (i. e, they are collections of lattice sites) on a 2D lattice, and interact with each other and the medium in a way that resembles cell-cell and cell-medium adhesion. Cells “move” by randomly copying their identity into neighboring lattice sites in order to minimize a global energy function, i. e. , the Hamiltonian H, that includes the cell-cell and cell-medium interactions and other cell constraints like cell area and cell perimeter: H = ∑ u, v J τ (σ u), τ (σ v) (1 - δ σ u, σ v) + ∑ σ λ Area (a σ - A σ) 2 + ∑ σ λ Perimeter (p σ - P σ) 2 where u and v are neighboring lattice sites, σ is the cell identity and τ is the cell type. The first term of the Hamiltonian sums the energy values, Jτ (σu), τ (σv), between all the neighboring lattice sites u and v, where the term 1 − δσu, σv (with Kronecker delta δσu, σv = 1 when σu = σv; and 0 otherwise) ensures that only interactions between different cells are considered in the Hamiltonian. The second term (the area constraint) and the third term (the perimeter constraint) sum over all cells; the system is penalized when the current cell area (aσ) deviates from the target cell area (Aσ) or when the current cell perimeter (pσ) deviates from the target cell perimeter (Pσ). The parameters λArea and λPerimeter are generally interpreted as resistance to compression and stiffness of the membrane, respectively. The area of a cell is the number of lattice sites with the same identity σ, while the perimeter is calculated as the number of distinct interfaces (edges and corners of lattice sites) with neighboring lattice sites of different cells or of the medium [39]. The lattice is updated asynchronously and one Monte Carlo step (MCS) is defined as the number of random copy attempts equal to the number of lattice sites. During each copy attempt, the algorithm picks a random site u and one of its neighbors v and, if they do not belong to the same cell (i. e. , σu ≠ σv), it tries to copy the cell identity of the site u into site v, effectively trying to allocate site v to the cell σu. In practice, the global value of H need not be calculated for every copy attempt; instead, the energy difference, Δ H, between the current configuration and the potential novel configuration is calculated from local information around u and v. A copy attempt is always accepted if it decreases the energy of the system (Δ H < 0), otherwise it is accepted with the Boltzmann probability e - Δ H / T, where T is the “temperature” of the system. Within the CPM, chemotaxis is implemented by increasing the probability of copying to sites with relatively high chemoattractant concentrations: Δ H Chemotaxis (u → v) = λ Chemotaxis (C v - C u), where Cu and Cv are the chemoattractant concentrations at lattice sites u and v. In order to contribute to the energy value of the system, Δ H Chemotaxis is subtracted from Δ H. In general, CPM can incorporate various other processes by subtracting additional terms from the Δ H. For an in-depth overview of CPM and its applications see [37]. We have implemented our model in the Tissue Simulation Toolkit, which is a two dimensional open source library for the CPM [19]. Additionally, an on-line interactive implementation is available at http: //bioinformatics. bio. uu. nl/ioana/cpm/ (tested in Firefox and Chrome browsers). All simulations are two dimensional and use a rectangular lattice. The CPM parameters that have constant values across simulations are the temperature T = 20, the area constraint λArea = 50 and the energy value of the CPM medium with itself, Jmedium, medium = 0. The perimeter constraint of all cells is λPerimeter = 2, with the exception of the tissue rigidity experiments where this value is varied for the tissue cells. In the multicellular simulations, we have used the connectivity constraint described by Merks et al. [19] to keep the Act cells from breaking. The migrating cells in the single cell simulations (Figs 2,3 and 4) and collective migration simulations (Fig 9) have a constant target area of Acell = 500, with λArea = 50 and a target perimeter of P = 340. The energy between cells is Jcell, cell = 100 and the energy between the cells and the medium is Jcell, medium = 20. The simulations without chemotaxis use a wrapped lattice of size 200×200; the simulations with chemotaxis use a wrapped lattice of size 100×300 and a linear chemoattractant concentration with a slope of 0. 33 applied in the y direction (i. e. , concentrations Cx, y = 0 = 0 and Cx, y = 300 = 100). In the collective migration experiments (Fig 9), we place Act cells of the same type on a wrapped lattice of size 500x500 and performed simulations in which we varied the cell type, the fraction of lattice coverage, and the cell-cell adhesion. To vary the cell type, we varied the parameter MaxAct, such that it covers the amoeboid to keratocyte-like range. We varied the fraction of lattice coverage between zero and one, where zero means that the lattice is empty and one that the lattice is full. To vary cell-cell adhesion, we varied the surface tension γcell, medium = Jcell, medium − (Jcell, cell + Jmedium, medium) /2; a positive γcell, medium encourages cells to “stick” together, while at γcell, medium = 0 cells adhere as much to the medium as to each other, i. e, a neutral value of adhesion. We varied γcell, medium by changing only the energy between the Act cells, Jcell, cell, while keeping Jcell, medium = 20 and Jmedium, medium = 0 constant. Jcell, cell = 40 results in γcell, medium = 0, the neutral adhesion value, Jcell, cell = −60 results in γcell, medium = 50, which we define as a medium adhesion, and Jcell, cell = −160 results in γcell, medium = 100, which we refer to as high adhesion. The target area of the Act cells is Areacell = 200 and their target perimeter λPerimeter = 180. In the skin model simulation (Fig 10), we created T cells by using the Act model parameters for amoeboid cells (MaxAct = 20, λAct = 2000). The target area of T cells is ATcell = 100 and the target perimeter is PTcell = 140. We created almost inflexible and almost stationary skin cells (the gray cells) by giving basic CPM cells (without the Act mechanism) a tight target perimeter (Pskin = 145) compared to their target area (Askin = 152), which prevents them from ruffling. The energy values in the simulation are Jskin, skin = Jskin, Tcell = Jskin, medium = Jskin, medium = 20, JTcell, Tcell = 100. For the exploratory tissue experiments (Fig 11), we performed simulations by placing one Act cell (MaxAct = 20 or MaxAct = 80 with λAct = 2000) in tissues with different densities, or in tissues with different rigidities at the highest density (a fraction of lattice coverage of 1). The target area and the target perimeter of tissue cells are the same as those of the skin cells in Fig 10, and the target area and target perimeter of the Act cells are the same as those of the T cells in Fig 10. The energy values in the simulations are Jtissue, tissue = Jtissue, ActCell = Jtissue, medium = JActCell, medium = 20, JActCell, ActCell = 100. We have varied the density of the tissue (Fig 11 left column) from relatively scattered cells (a fraction of lattice coverage of 0. 25) to a full lattice (a fraction of lattice coverage of 1). We have varied the rigidity of the tissue (Fig 11 right column) by varying the value of λPerimeter of the tissue cells. This parameter governs how stringent a cell conserves its target perimeter: at λPerimeter = 0 a cell ignores the perimeter constraint, while at a higher value of λPerimeter the cell is allowed less deviations from the target perimeter. A higher value of λPerimeter combined with a relative short target perimeter creates cells that are more difficult to deform, hence they are more rigid. Centroids and principal axes of cells are calculated from the moments of inertia [40]. The length of a cell is the length of its major axis. The instantaneous speed of a cell is the Euclidean distance traveled by its centroid during one MCS. The orientation-direction angle (∈ [0. 90]) is the angle between the cell’s longest axis and its direction of migration. The turning angle is the angle between two consecutively measured directions of migration. Persistence time is expressed as the number of MCSs required for a cell to lose the angular correlation of the direction of migration [41]. The motility coefficient is calculated by fitting Fürth’s equation for persistent random migration [42,43] to mean square displacement plots: x 2 ¯ = 2 n M (t - P (1 - e - t / P) ), where x 2 ¯ is the mean square displacement, n is the spatial dimension (here n = 2), M is the motility coefficient, P is the persistence time and t is the elapsed time since the start of the trajectory. For the analysis of the cell chemotaxis simulations, we use the directed and scaled directed speed as measurements for chemotaxis sensitivity. The directed speed is defined as the magnitude of the cell velocity projected onto the chemokine gradient. The chemotactic index is defined as the distance traveled up the chemokine gradient divided by the total distance traveled. A cell protrusion is defined as a connected patch of at least 10 active lattice sites. The orientation of a protrusion relative to the chemokine gradient is defined as the angle between the chemokine gradient and the vector that connects the centroid of the cell with the centroid of the protrusion. The collectivity of cell migration is calculated as the order index, | ∑ i = 1 n v i → | v i | |, where n is the number of cells on the lattice and v i → is the direction vector of cell i during one MCS. In order to ensure the evaluation of collectivity at steady state, the order index of one simulation is calculated at the end of 500. 000 MCSs.
Cell migration is involved in vital processes like morphogenesis, regeneration and immune system responses, but can also play a central role in pathological processes like metastasization. Computational models have been successfully employed to explain how single cells migrate, and to study how diverse cell-cell interactions contribute to tissue level behavior. However, there are few models that implement realistic cell shapes in multicellular simulations. The method we present here is able to reproduce two different types of motile cells—amoeboid and keratocyte-like cells. Amoeboid cells are highly motile and deform frequently; many cells can act amoeboid in certain circumstances e. g. , immune system cells, epithelial cells, individually migrating cancer cells. Keratocytes are (fish) epithelial cells which are famous for their ability to preserve their shape and direction when migrating individually; during wound healing, keratocytes migrate collectively, in sheets, to the site needing reepithelialization. Our method is computationally simple, improves the realism of multicellular simulations and can help assess the tissue level impact of specific cell shapes. For example, it can be employed to study the tissue scanning strategies of leukocytes, the circumstances in which cancer cells adopt amoeboid migration strategies, or the collective migration of keratocytes.
Crawling and Gliding: A Computational Model for Shape-Driven Cell Migration