ORIGINAL RESEARCH article
Modeling Active Cell Movement With the Potts Model
- 1INIFTA, Universidad Nacional de La Plata - CONICET, La Plata, Argentina
- 2Instituto de Investigaciones Científicas y Tecnológicas en Electrónica, Universidad Nacional de Mar del Plata - CONICET, Mar del Plata, Argentina
- 3Centro Regional de Estudios Genómicos, Universidad Nacional de La Plata - CONICET, La Plata, Argentina
In the last decade, the cellular Potts model has been extensively used to model interacting cell systems at the tissue-level. However, in early applications of this model, cell movement was taken as a consequence of membrane fluctuations due to cell-cell interactions, or as a response to an external chemotactic gradient. Recent findings have shown that eukaryotic cells can exhibit persistent displacements across scales larger than cell size, even in the absence of external signals. Persistent cell motion has been incorporated to the cellular Potts model by many authors in the context of collective motion, chemotaxis and morphogenesis. In this paper, we use the cellular Potts model in combination with a random field applied over each cell. This field promotes a uniform cell motion in a given direction during a certain time interval, after which the movement direction changes. The dynamics of the direction is coupled to a first order autoregressive process. We investigated statistical properties, such as the mean-squared displacement and spatio-temporal correlations, associated to these self-propelled in silico cells in different conditions. The proposed model emulates many properties observed in different experimental setups. By studying low and high density cultures, we find that cell-cell interactions decrease the effective persistent time.
Biological development is a clear example of a complex system and is the result of several morphogenetic processes such as: cell adhesion, cell division, migration, among others . Important insights about the principles involved in the biological development can be extracted from mathematical models used in many works over the last decade . However, there is still no standard model to describe interacting cells in all scales. A suitable approach depends on the spatial and temporal scale associated to the question under study [3, 4]. Continuous models to study combinations of chemical reaction kinetics and cell mobilities through diffusional processes have been useful at tissue-level scales [5, 6]. Also, simplified continuous models of self-propelled particles have been useful to uncover a new kind of phase transition in collective motions [7–12], (see  for a review of motility-induced phase separation).
However, these models could not be adequate to analyze biological processes involving adhesion and/or changes in cell shapes. On the other hand, cell-based models, such as the lattice-gas cellular automaton [14, 15] or cellular Potts models (CPM) [16, 17] offer a better choice to analyze interacting cell systems at the tissue-level. These discrete models usually do not incorporate physical-chemical processes in detail, since it would imply covering a broad range of scales, and are not amenable to mathematical analysis. However, particularly CPM has been shown to be suitable for studying very complex phenomena, via computer simulations. In this sense it has been applied successfully in the study of cell sorting in cellular aggregates [16, 18], morphological development of the Dictyostelium discoideum slug [19–21], and tumor growth and angiogenesis [22, 23]. Recently, a continuum limit of the CPM was proposed to study the role of membrane elasticity on the migration potential of cancer cells .
CPM was introduced by Glazier and Granier in '90 [16, 18], as an extended large-Q Potts model. CPM considers individual cells, which can suffer small deformations and volume changes according to a free energy minimization principle which drives its dynamics. As a consequence of the small changes resulting from cell interactions with elements of their immediate surrounding, the cells perform small displacements (≪ cell size ). Statistical properties of this “passive" movement of in silico cells in aggregates have also been studied [26, 27]. With this model, cells perform a random walk; i.e., the position and velocity power spectra correspond to Brownian motion and the velocities are time-uncorrelated. Although some correlations have been detected both in Hydra cell aggregates  and in Madin-Darby canine kidney cells [29, 30], this finding is in agreement with experiments performed with pigmented retinal cells .
Cell motion is usually influenced by processes that involve several components of the cellular migration machinery, resulting in cell polarization. When the direction of polarization is maintained in time the cell can perform a persistent movement. This active movement has been observed in collective migration both in response to external chemotactic gradients , as well as in the absence of external signals [29, 30, 32–34]. Persistent polarization, intermittency and changes in the movement direction are consequently reflected in strong spatial and temporal correlations [29, 30]. In order to direct cell motion in the CPM framework, researchers have added terms to the differential of the Hamiltonian that drive a directional bias of cell movement. Thus, a term has been introduced to model the chemotactic response of Dictyostelium cells in response to cAMP signaling during the slug morphogenesis . Similar bias terms have been proposed for making cell motility dependent on a propensity vector, thus modeling persistent, active cell motility in absence of external fields [33, 36].
In the same direction, it has been used more recently in the context of collective cell migration or active matter [37, 38]. These simplified models consider persistent cell motion externally embedded in the CPM. In contrast, some efforts have been devoted to build models that include an actin-based mechanism resulting in cell protrusions and, consequently, in active cell movement , or even more detailed models that integrate the biochemical reactions and cytoskeleton dynamics to cell motion . These multiscale models can be used to understand how molecular-based mechanisms give rise to the epiphenomenon of persistent cell motion in the absence of external signals.
In this paper, we study the effects of low and high density cultures in persistent cell movement when external signals are absent. For this purpose we use the CPM, which considers cell adhesion and deformation, to better describe cell-cell interactions. Alternatively, one can use a particle-based model with volume exclusion, but in our opinion this is not the best way to consider interaction between cells. Therefore, we propose a modification in the dynamics of the directional propensity vector introduced in previous works that consider active cell movement in the absence of external cues [33, 36]. In this sense, we link the polarization direction of the cell to an autoregressive process that changes the migration direction after some characteristic time τ, and study the dynamics of the resulting cell movements. We perform computer simulations for different parameter values and find that the dynamics of these in silico-cell movements resembles the behavior observed experimentally in the absence of external signals [34, 41]. We believe that the present modeling approach can fill the gap between interesting characteristics of the cell motility (correlations, velocity distributions) and key ingredients, such as persistence time and cell density.
2.1. The Model
The CPM has been introduced as an extension of the Potts model, which is a generalization of the Ising model. The CPM has been used in literature by many authors to study a wide range of phenomena, ranging from critical behaviors of bubbles to cellular arrangements . The model basically consists of a system of interacting spins σ on a lattice, where the set of all (adjacent) sites with the same spin defines a cell. Different terms of energy can be considered according to the purpose of each case study. In the basic model, the energy E0 is associated with the cell-cell interaction, cell volumes , and cell perimeters . Fluctuations of the membrane allow the cells to explore their neighborhood. Thus, if the total number of cells is Q, all sites i where σi = M, with 1 ≤ M ≤ Q, belong to the cell labeled as M. The energy of any configuration of cells is given by
where the primed sum runs over neighboring site pairs and δσiσj is the Kronecker delta, and equals one if σi = σj and zero otherwise. Thus, the first sum is over all neighboring site pairs of adjacent cells, and represents the boundary energy of the interacting cells. The second and third sums represent the energies of area and perimeter elasticities of the cells, respectively, and are related with the constants κ and Γ. In this way, a target area A0 and a target perimeter L0 are introduced and the sums are over all cells on the lattice, being AM and LM the area and perimeter of the cell M. The second and third terms on Equation (1) are required to generate domains of spins resembling cells, that will not break nor disappear, in a sustained manner over time . As we will see later, the presence of a medium, which interacts with the cells, can also be considered by assigning zero value to the spin variable σi = 0. Obviously, the medium has neither target area nor target perimeter.
The system evolves using standard Metropolis algorithm, with the following variant: a lattice site i is randomly chosen and if it belongs to the boundary of the cell M, then the site assumes the spin value of its neighboring cell M′ with probability proportional to the Boltzmann factor , where ΔE0 is the change of energy involved in the proposed update and kBT simulates the energy associated to membrane fluctuations driven by the cytoskeleton . At each Monte Carlo step, the total number of sites is randomly chosen and the probability of the configuration update is evaluated in each case. Thus, according to Equation (1), changes in cell shape and cell displacement are a result of membrane fluctuations alone.
The energy above has been broadly used in CPM to drive the cell movement [26, 42–44]. However, the terms in Equation (1) are insufficient to describe active cell motility, characterized by a persistent direction of cell movement. In order to simulate the active cell movement we consider an additional term in the energy, promoting a bias in the cell extensions toward a given direction, as originally proposed in Savill and Hogeweg . This term has been broadly used in modeling active cell motion both in the presence and absence of external cues [21, 33, 36, 37, 45], acting as a source of energy that drives motility of a given cell M preferentially along the direction of a driving force . Thus, the overall energy variation at each Monte Carlo step can be written as:
where ΔE0 is the change of energy due to Equation (1) and represents the displacement of the center of cell M.
The driving force is characterized by a strength and a direction denoted by ΘM. Here, while the strength is equal for all cells and constant over time, the driving force direction is different for each cell and evolves independently of the others. In fact, ΘM is actualized according to the following first-order autoregressive process:
where t* refers to the time of actualization of ΘM, which is different from the Monte Carlo time step used to update the configuration of the cells. ϕ and ΔΘ are constant parameters, that can take values between (0, 1) and (0, π), respectively. ϵ(t) is a white noise with zero mean and unit variance (). The stochastic difference Equation (3) can be reinterpreted in terms of the Ornstein-Uhlenbeck process with media zero, dx(t) = −λx(t)dt−σdW(t), where dW(t) denotes the Wiener process. In fact, setting the friction coefficient λ = 1−ϕ and σ = ΔΘ, Equation (3) becomes the discrete-time counterpart of the Ornstein-Uhlenbeck process. Over time, the mean value and the variance of Θ are 〈Θ〉 = 0 and , for ϕ < 1, respectively. On the other hand, the temporal autocorrelation of the angle ΘM, defined by the quotient between the autocovariance and the variance σΘ, is given by
Thus, the expected overall time decay for a cell movement in a given direction is −τ/ln ϕ. When ϕ closes to 0, the time decay tends to zero and the angle dynamics looks like white noise. But when ϕ approaches 1, and the temporal autocorrelation is large. Equation (3) expresses a putative dynamics for the driving force direction. Different update rules have been considered in previous works [33, 36, 37]. In fact we also explore an alternative one, that includes a positive feedback loop, in the Results section.
The movements of the cells do not necessarily have the same direction of the driving force, due both to the stochastic fluctuations and to cell-cell interactions. For example, if the driving force is applied on two adjacent cells in opposite directions and one of them attempts to move in the direction of the other, there is no contribution to the energy variation. In order to avoid confusion, we refer to the angle of the driving force and of displacement of the cell M as ΘM and αM, respectively.
2.2. Simulation Details
To compute the energy contribution associated with cell-cell and cell-medium interactions we consider a neighbor range greater than one lattice site. Thus, the primed sum in Equation (1) considers 1st, 2nd, 3rd, and 5th neighbors, which contribute to the energy Jσiσj with weights proportional to the inverse of their distances di,j. In short, if σi ≠ σj ≠ 0, the surface energy Jσiσj = Jcell−cell/di,j, where Jcell−cell is the cell-cell adhesion constant. On the other hand, if one of the neighboring sites, i or j, belongs to the medium, the surface energy is given by Jσiσj = Jcell−medium/di,j, where Jcell−medium is cell-medium adhesion constant. Therefore, the first term of Equation (1) contributes to the energy if the sites i and j belong to different cells or one of them belongs to a cell and the other to the medium, representing the cell-cell or cell-medium adhesion, respectively. As cells tend to stick together, we chose Jcell−cell > Jcell−medium in our simulations, as we will see next.
The initial configuration is generated by partitioning the lattice into equal-sized square domains, corresponding to the cells, each cell containing 16 × 16 lattice sites. The squares alternate offsets in every other row, so the pattern resembles a brick wall arranged in common bond. This situation corresponds to a dense aggregate of cells which have been adopted in most previous studies. As the diameter of the cell is about 16 pixels, the parameter L0 and A0 have been set to π × 16 and π × 64, respectively. For comparison purposes, taking into account that most eukaryotic animal cells have a diameter between 10 and 30 μm, we can establish the equivalence μm = pixel.
In order to understand the role of cell-cell interactions on the statistical characterization of locomotion activity, we also study cellular movement at low and high density cell cultures. To this end, we assign the same spin value to a number of randomly chosen cells, σi = 0 associated with the medium. Thus, the zero spin domain represents the medium, whereas the group of spins σi = M with 1 ≤ M ≤ Q values, belong to the cell labeled as M. The density ρ is defined as the quotient between the area occupied by cells with σi ≠ 0 and the total area, thus ρ = 1 implies a lattice fully populated with cells.
Simulations start with a given cell density ρ, allowing for a “thermalization” time interval on the dynamics to obtain adequate cellular shapes, during which no motility forces are applied on cells. In this thermalization interval, of the order of 100 MCS, the configuration is updated through calculating the Boltzmann factor according to Equation (1). Once the system thermalizes, we establish the starting time t = 0 (see Figure 1). At this point the dynamics starts, with the presence of the motility force on the cells with constant strength and angle given by Equation (3). The initial angle of each cell, ΘM(0), is chosen randomly. At each Monte Carlo step the driving force direction for each cell changes with probability 1/τ according to Equation (3), where indicates the direction previous to the actualization. Thus, the change in the driving force direction ΘM occurs at a mean time τ.
Figure 1. Typical snapshots of a simulation carried out on a lattice of 512 × 512 sites at two different times: (A) initial time (t = 0) in which the dynamic starts with the presence of driving forces acting on each cell, (B) after t = 2, 000 MCS the cells have moved on the medium. For example, gray colored cell has displaced across scales much larger than the cell size. The parameter values used for this simulation were ρ = 0.2, τ = 10 MCS, and ϕ = 0.95.
For all simulations in this paper, we used periodic boundary conditions and a square lattice of size 512 × 512 sites. We have also fixed the parameter values Jcell−cell = 0.1, Jcell−medium = 0.01, Γ = 0.2, κ = 1, ΔΘ = π/3, F = 10, and T = 2; while we vary the values of parameters ρ, τ and ϕ in order to study the dynamics of the resulting cell movement. Figure 1 illustrates two instances of the simulation at low density culture (ρ = 0.2), with a driving force of constant module and whose direction is governed by Equation (3). Additionally, the dynamics of the resulting simulation for low and high density cultures can be seen in the Supplementary Movies 1, 2, respectively.
Our selection for ΔΘ avoids a preferential rotation of cells, since the stochastic term in Equation (3) blurs the angle of the previous time. Then, cells perform random walks in the long time scale as we will show in subsection III-C. For small values of ΔΘ cells can display a preferential rotation, which has been observed for cells confined in circular micropatterns [38, 46]. However, the study of geometrical confinement into well-defined circles is outside the scope of this work.
Figure 2 shows four different simulated cell trajectories and their associated angles on a 2D substrate using our model with the same parameter values as those used in Figure 1 A combination of persistent and random cell displacements can be seen in Figure 2A. These trajectories resemble those traced by Dictyostelium cells , which travel much longer distances than the typical cell size. It is important to note that a constant driving force, when applied to an isolated cell, originates a movement of constant velocity . In this way, the persistent movement in our model is due to the driving force, that has a constant strength and is applied with constant direction during a mean time τ, as discussed before.
Figure 2. (A) The trajectories are marked off with red circles each 100 MCS and the squares correspond to the positions of each cell at t = 0 MCS. (B) Temporal course of the cell displacement angles α associated with the trajectories depicted on the left panels. Parameters used: ρ = 0.2, τ = 10 MCS, and ϕ = 0.95.
The panels of Figure 2B depict the cumulative cell displacement angle and were obtained by tracking the cell displacement angles with reference to a fixed system. In general, over a short time scale, cells move in a certain direction around which the instantaneous angle direction fluctuates. These fluctuations are due to the cell membrane movement driven by the cytoskeleton and also by cell interactions. At large time scales (>τ) random changes in the cell direction become more apparent and they are promoted by the driving force whose angle actualization is given by Equation (3). Over time, the cumulative cell displacement angle, in addition to offering the instantaneous direction of the cell at a given time, also tracks the entire history of the cell's directional changes. Therefore, periods of persistent cell displacement correspond to constant average angles, while large turns in the cell trajectories are related to large increases or decreases in the average angles. In this way, cell 4 maintains a relatively persistent movement between 30 and 200 MCS (Figure 2A) which is accompanied by a cell displacement angle that fluctuates around the constant value α~−0.8πrad (Figure 2B); and cell 2 presents marked increases and decreases in α around the first 100 MCS (Figure 2B) that correspond to different turns in its trajectory (Figure 2A).
3.1. Cell Displacement
We studied several statistical properties associated to the cell displacement which can be compared with the ones derived from experimental measures. In this sense, we define the cell displacement of a cell M, in a time interval Δt, as . Associated with this cell displacement, there is a displacement angle αM(t) and a mean velocity during this interval.
First of all we consider the cell displacement angle α by using Δt = 1. As shown in the histogram of Figure 3A, cells can move in any direction with similar probabilities, putting in evidence that there is not a preferential direction. This histogram was computed over all cells during a simulation with the same parameter values as those used in 1. Figure 3B depicts the temporal autocorrelation function of α, defined as Cα(t) = Zα(t)/Zα(0) where Zα(t) = 〈αM(t0 + t)αM(t0)〉, and the averaging was performed over all cells. Besides, the dynamics of α has a characteristic time, which depends on the parameters τ and ϕ, as suggested by expression (4) for Θ. The characteristic time computed in the exponential decay region (inset of Figure 3B) is 194 MCS, in agreement with the theoretical prediction associated to the driving force angle, −τ/ln ϕ = 195 MCS.
Figure 3. (A) Frequency distribution of cell displacement angle α for all cells during 4,000 MCS. (B) The autocorrelation of α, Cα, as a function of time. Inset: log scale of Cα showing the region of exponential decay (short times). (C) 2D-histogram of the turn angle (Δαt versus Δαt + 1), showing that two consecutive turn angles are anti-correlated. Parameters used: ρ = 0.2, τ = 10 MCS, and ϕ = 0.95.
Figure 3C depicts a 2D-histogram of the recurrent plot of the turn angles Δθ(t) vs. Δθ(t + 1), where Δθ(t) = θ(t)−θ(t−1). This plot shows that two consecutive turn angles are anti-correlated. This feature of cell motility has also been observed in experiments with Dictyostelium cells .
As observed by Cox , some statistical properties associated with cell movements depend on the time interval Δt considered. We studied the distribution of the cell velocity components vx and vy for different time intervals Δt = 1, 3, 10, and 60 MCS. Figure 4 illustrates the distribution of cell velocities for ϕ = 0.95, τ = 10 and two different densities: ρ = 0.2 (Figure 4A) and ρ = 0.9 (Figure 4B). These distributions were obtained by tracking 51 cells (chosen at random). For Δt = 1, the velocity distributions for both densities are quite similar, resembling Gaussian distributions, but with small gaps in the centers of the plots. However, as Δt increases, the difference between the distributions of the cell velocity components obtained at low and high densities becomes more apparent. When Δt approaches to τ, the gap in the center increases for ρ = 0.2, but not for ρ = 0.9, where it is quite smaller and only apparent for Δt = 3. In fact, the difference between low and high densities is more evident when Δt = τ, indicating that this is an adequate time interval size to characterize the underlying dynamics. Similar results were obtained from simulations by using τ = 50 (Supplementary Figure 1), where the main difference between high and low densities is also found for Δt = τ.
Figure 4. Cell velocities vx vs vy calculated for different values of the time interval Δt (upper panels). Histogram of the cell velocity component vx computed over the vy window indicated by red bars (lower panels). Two different values of densities were used: ρ = 0.2 (A) and ρ = 0.9 (B). All other parameters are equal to Figure 2.
It is important to note that the gap in the velocity distribution was not observed in simulations performed when the persistence time is absent (i.e., τ = 1, Supplementary Figure 2). Thus, the presence of a gap in the center of the velocity distribution is associated with effective persistent cell movement, in agreement with Maiuri et al. (Figure 4D of ). On the other hand, in high density simulations the crater-like effect is less evident (Figure 4B and Supplementary Figure 1). Note that for ρ = 0.9 and Δt = τ the velocity distribution resembles a Gaussian one, indicating that cell-cell interactions, massively present in these cases, contribute to decrease the effective persistent movement of the cells. Large enough values of Δt (>τ) blur the instantaneous direction of cell, and distributions are characterized by a bell shape independently of the densities. Similar results were obtained for ϕ = 0.99 (Supplementary Figure 3), where the behavior of the distributions tails for ρ = 0.9 are more Gaussian than for ϕ = 0.95.
The movement of Dictyostelium cells in the absence of external signals presents a sequence of velocity distributions by varying Δt  similar to that obtained from our model, for ρ = 0.2.
3.2. Velocity Correlation
Many studies on cell motility have considered both temporal and spatial correlation functions of the cell velocity for statistical characterization of locomotion activity [25, 28–30, 34]. The temporal autocorrelation function is defined as C(t) = Z(t)/Z(0) where . The cell velocities were computed using Δt = 1 MCS, and 〈…〉 indicates the average over all cells and over t0.
Figure 5A depicts the behavior of C(t) at low density (ρ = 0.2), for three values of the parameter τ = 10, 50, and 100, and two values of ϕ = 0.95 and 0.99. As expected, the correlation decay increases with τ. For ϕ = 0.99 (dashed lines) C(t) exhibits a single exponential decay, which is consistent with an Ornstein-Uhlenbeck process. However, for ϕ = 0.95 (solid lines) C(t) exhibits a more complex behavior: a double exponential decay with a slower decrease at long time scales, independent of the τ value. Moreover, as expected, the autocorrelation function associated to the driving force with ϕ = 0.95 is greater than that obtained from the more random driving force, ϕ = 0.99. We also notice that the behavior of C(t), obtained with the same ϕ but different τ values, coalesces when it is scaled with the parameter τ (Figure 5C).
Figure 5. Semi-log plot of the temporal autocorrelation function C(t) averaged over all cells in the simulation. The plots correspond to three different values of τ, two values of ϕ (0.95 solid lines and 0.99 dashed lines) and for two different densities: ρ = 0.2 (A) and ρ = 0.9 (B). (C,D) Shows the collapse of the curves by rescaling the horizontal axes by τ, for the simulations at low and high density, respectively. All other parameters values are the same as in Figure 2.
However, the effects of the cell-cell contact on the autocorrelation function are much more complex. Similar simulations to those used in Figure 5A but at a higher density, ρ = 0.9, show that C(t) does not exhibit a single exponential decay in none of the cases studied (Figure 5B). On the other hand, the effect of the density on the decay of the autocorrelation function depends on the value of the persistence time τ. In this sense, when comparing short and long-time behavior, we can observe that the decay rate increases at long time scales for τ = 10, while it decreases for τ = 100. Consequently, the behavior of C(t) for different τ values does not coalesce when scaled with τ (Figure 5D). We found that the temporal autocorrelation obtained at high density cultures is smaller than at low density cultures. This fact can be understood in terms of the decrease of the effective persistence time of cell movements due to a crowded neighborhood, as seen also in Figure 4.
Spatial correlations of the cell velocities have also been used to characterize the dynamics [25, 28, 29]. The spatial correlation is defined as , with the distance between mass center of cells M and M′ . Figure 6 shows C(r) for two different cell densities on the substrate, ρ = 0.2 and 0.9 (solid and dashed lines, respectively), and three different values of τ. These curves display a peak around 14–16 pixels, close to the diameter of a non-interacting cell. Since the in silico cells have an elastic membrane, they can be deformed to reach intercellular distances smaller than the cell diameter. This can occur when the cells travel in opposite directions or in order to increase the contact area. In Figure 6 one can observe that at very short distances (r < 8 pixels) cells are anti-correlated. The negative values for the spatial correlation correspond to cell-pairs in collision, i.e., cells traveling in opposite directions. These cell-pairs are also present, in less extent, in the range 8 ≤ r ≲ 14 pixels, contributing to lower the correlation. Thus, in this range, one can observe that C(r) obtained for simulations at low density are higher than the ones obtained at high densities, in particular for τ = 10. However, for longer distances (r ≳ 16 = cell diameter) this relationship is inverted, i.e., cells in high density configurations are more correlated. Besides, C(r) in high density simulations approaches zero very slowly, as r increases. These results are independent of the τ-value, indicating that the cell-cell contact induces long-range spatial-correlation of cell velocity. Interestingly, highly correlated motions have been found in crowded tissues  and in cells in different aggregate types . Experimentally, high spatial correlations have been found in cell aggregates moving together collectively .
Figure 6. Average spatial correlation function of the velocities, C(r), as a function of the distance r between cell pairs, for two different densities, ρ = 0.2 and 0.9 (continuous and dashed lines, respectively) and different values of τ, as indicated, and ϕ = 0.95. The vertical dashed line represents the average cell diameter equal to 16 μm as discussed in section 2.2. Data was obtained by averaging between 0 and 4,000 MCS and over 50 and 10 samples, for ρ = 0.2 and 0.9, respectively.
At intermediate distances C(r) decays exponentially allowing to define a correlation length. We found that, the correlation lengths obtained for ρ = 0.9 are larger than those obtained using ρ = 0.2, independently of the τ-value. Also, for ρ = 0.9, simulations with τ = 10 present a slightly smaller correlation length than simulations with τ = 50 and 100 (9.1, 10.2, and 11.0, respectively). This means that correlation length increases with cell velocity, as cell velocity also increases with τ. This fact is in agreement with previous findings reporting the increase of correlation length with cell velocities . On the other hand, for ρ = 0.2, all simulations present similar values for the correlation length (5.2, 5.4, and 5.4 for τ = 10, 50, and 100, respectively).
Note that, according to Equation (3) the angle of the driving force ΘM depends on its previous value but does not depend on the direction of the cell movement αM. Other authors have considered the presence of a positive feedback loop in cell polarization dynamics, in which cell produces its own chemotactic signal [33, 36, 37]. In order to take into account this feedback loop, an alternative way to update the angle ΘM is:
where is the mean displacement angle over the last τ MCS, and the other parameters are the same as in Equation (3). Thus, according to Equation (5) the angle of the driving force depends on the earlier cell displacements with a weight given by the parameter ϕ.
We observe that, for any cell density, both the temporal and spatial correlations of the cell velocities are higher when using Equation (5) instead of Equation (3) (see Supplementary Figures 4, 5). In particular, the increment in the correlations is more evident for ρ = 0.9, indicating that the differences between the angle update rules are more apparent in crowded cell cultures. C(t) obtained when considering Equation (5) for ρ = 0.2 is almost independent of the parameter ϕ, whereas for ρ = 0.9 short and long-time behavior of C(t) depends on τ and ϕ in a non trivial manner (see Supplementary Figure 4). Again, the correlation length obtained with Equation (5) increases with τ, reaching 1.6-fold of the cell size when τ = 100 (see Supplementary Figure 5). However, this value is smaller than the correlation length of 2–3 cell size, reported in Garcia et al. . Interestingly, when a feedback from previous displacements is taken into account through Equation (5), C(r) is less negative at very short distances, suggesting that adjacent cells can be avoided in a more realistic way.
3.3. Mean-Squared Displacement
The mean-squared displacement (MSD) is a useful measure to find cell persistence of migration in two and three dimension microenvironments [49–52]. This is calculated as , where 〈…〉 indicates the average taken over all the cells on the substrate. Figure 7 shows MSD scaled with time as a function of time at ρ = 0.2 (solid lines) and at ρ = 0.9 (dashed lines) and three different values of τ. The MSD obtained for high density is lower than that obtained for low density, as expected. In all cases, MSD behaves as ~t2 for short enough time scales (ballistic), and MSD behaves as ~t for long time scales (random walk), as expected for cell movements in the absence of symmetry-breaking gradients [34, 41]. Scaling MSD with τ2 and the time with τ, we can observe that the scaled behavior is independent of ϕ and τ for ρ = 0.2. However, for high density, there is a small deviation and scaled behavior depends on the τ-value, showing the complex behavior generated by the interaction between cells.
Figure 7. (A) Log-log plot of the scaled mean-squared displacement vs. time t, for three different values of τ, ϕ = 0.95 and for two different densities: ρ = 0.2 (solid lines) and ρ = 0.9 (dashed lines). Vertical arrows mark the transition between ballistic and random walk behavior. (B,C) shows the collapse of the curves by rescaling the horizontal axes by τ and vertical axes by τ2, for the simulations at low and high density, respectively. All other parameters values are the same as for Figure 2.
Since the classical Steinberg experiments, it has been demonstrated that differential adhesion is a key ingredient for cell rearrangement and cell sorting [53, 54]. Graner and Glazier have shown that the CPM model with solely differential adhesion, as driving force, is able to generate cell sorting . However, during the morphogenetic processes, there are other ingredients that also play important roles, like chemotaxis. Savill and Hogeweg have extended the CPM by considering a “directional propensity" vector to model the chemotactic response of Dictyostelium cells . This vector can induce a bias on the membrane sites belonging to a cell in the direction given by chemoattractant gradients. Thus, the cells move at a constant velocity whose module depends almost linearly on the intensity of the bias . Accordingly, the extension of the CPM allows to leave considering cells as passive objects, giving them the possibility to actively explore the environment beyond the local fluctuations of their cell membranes, both in the presence or absence of external signals. This extension has been used to model the D. discoideum morphogenesis [19–21], to investigate collective cell migration (or motion) in Szabó et al. , Kabla , Czirók et al. , and Czirók et al. , and more recently to study the chemotactic response of cells . In this paper, following the line of pioneering works that use the CPM to study persistent cell motion in the absence of external cues [33, 36], we take advantage of a directional propensity vector linked to a first order autoregressive process that governs the dynamics of cell movement direction. Thus, the timing of the cell movement is given by two parameters: the persistence time τ and the autoregressive parameter ϕ, instead of only one as in Beltman et al.  and Szabó et al. .
We study here some statistical features of the cell movement in this model under different conditions and observe that it is able to reproduce several features of the motility of real cells observed in the absence of chemoattractants. It has been pointed-out that the motion of Dictyostelium and Polyphondylum cells does not follow the Brownian motion and, unlike the Ornstein-Uhlenbeck process, cell velocity distribution deviates from Gaussian distribution . In the same line, we observe a gap in the velocity distribution when an effective persistence time is present. On the other hand, bell-shape distributions were observed when τ = 1 or in a crowded neighborhood. Both bell- and crater-shape distributions for cell velocity have been predicted by a mathematical model that considers the actin flow , connecting persistent motion with crater-shape distribution, in agreement with our findings. Also, we observe in our simulations at low density cultures that the temporal autocorrelation of the cell velocity obtained for ϕ = 0.99 decays exponentially, while for high density and/or ϕ = 0.95 the behavior deviates from the expected random motion under the traditional Ornstein-Uhlenbeck process. For ϕ = 0.95 and low density, we observe a double exponential decay with a slower decrease of the autocorrelation function at long time scales independent on the τ value. This two-step profile for the autocorrelation function has been previously observed in different experimental setups [41, 57]. On the other hand, the effect of the density on the autocorrelation function depends on the value of the persistence time. Moreover, similarly to the effect of ϕ = 0.95, the cell-cell interaction induces a no-single exponential decay behavior. This complex decay behavior could be the result of the interplay between the parameters τ and ϕ. Besides, we find that the particular way in which the cell movement direction is updated has important effects on the spatial correlation. In this sense, when this update rule considers a feedback between polarity and cell motion, similarly to proposed in Szabó et al. , we obtain larger correlation length than when this feedback is absent. We also find that high density cultures have larger correlation length than low density cultures, suggesting that cell-cell interactions play a crucial role in cell collective motion. Summing up, these results suggest that the representation of individual cells and their interactions, allowed by the CPM, can be critical to study the effect of different density cultures on the statistical feature of the cell movement.
While the parameter ϕ has important effects over the temporal autocorrelation function, almost no effect can be seen in the mean-squared displacement, mainly modified by the cell-cell interactions. We found that these results are consistent with the prediction of the persistent random walk model [49–51]. This agreement is not in itself enough to completely specify the dynamics of cell motion, since other processes used to model random motion also share this feature [25, 58]. The behavior of the mean-squared displacement observed in our simulations is in agreement with the one observed in different experimental setups [34, 41]. Moreover, the temporal autocorrelation function of the cell velocities also deviates from the single exponential decay predicted by the persistent random walk model. It is interesting to note that anomalous diffusion of cells has also been reported from different experiments [29, 30]. This type of cell motion is better understood in terms of the fractional Klein-Kramers equation, which includes the Ornstein-Uhlenbeck process as a special case .
These examples put in evidence the difficulty to find in a single model the ability to accurately represent the full range of cell movements observed in nature. Therefore, the model proposed here can be added to the cell movement model repository with promising perspectives.
NG, KM, and LD conceived the paper, revised the manuscript and approved the final version. LD developed the the structure and design of the manuscript. NG and KM developed the algorithms and wrote the manuscript.
This research was supported by CONICET, Argentine (PIP No. 0597).
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
We would like to thank Luciana Bruno and Fernando Peruani for the fruitful discussions and for their helpful comments on this work.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2018.00061/full#supplementary-material
9. Sepúlveda N, Petitjean L, Cochet O, Grasland-Mongrain E, Silberzan P, Hakim V. Collective cell motion in an epithelial sheet can be quantitatively described by a stochastic interacting particle model. PLoS Comput Biol. (2013) 9:e1002944. doi: 10.1371/journal.pcbi.1002944
11. Peruani F, Starruß J, Jakovljevic V, Søgaard-Andersen L, Deutsch A, Bär M. Collective motion and nonequilibrium cluster formation in colonies of gliding bacteria. Phys Rev Lett. (2012) 108:098102. doi: 10.1103/PhysRevLett.108.098102
17. Alber MS, Kiskowski MA, Glazier JA, Jiang Y. On cellular automaton approaches to modeling biological cells. In: Mathematical Systems Theory in Biology, Communications, Computation, and Finance, New York, NY: Springer (2003), 1–39.
20. Marée AFM, Hogeweg P. How amoeboids self-organize into a fruiting body: Multicellular coordination in Dictyostelium discoideum. Proc Natl Acad Sci USA. (2001) 98:3879–83. doi: 10.1073/pnas.061535198
27. Ouchi NB, Glazier JA, Rieu JP, Upadhyaya A, Sawada Y. Improving the realism of the cellular Potts model in simulations of biological cells. Phys A (2003) 329:451–8. doi: 10.1016/S0378-4371(03)00574-0
28. Upadhyaya A, Rieu JP, Glazier JA, Sawada Y. Anomalous diffusion and non-Gaussian velocity distribution of Hydra cells in cellular aggregates. Phys A (2001) 293:549–58. doi: 10.1016/S0378-4371(01)00009-7
29. Diambra L, Cintra LC, Chen Q, Schubert D, Costa L da F. Cell adhesion protein decreases cell motion: statistical characterization of locomotion activity. Phys A (2006) 365:481–90. doi: 10.1016/j.physa.2005.10.006
31. McCann CP, Kriebel PW, Parent CA, Losert W. Cell speed, persistence and information transmission during signal relay and collective migration. J Cell Sci. (2010) 123:1724–31. doi: 10.1242/jcs.060137
40. Marée AFM, Jilkine A, Dawes A, Grieneisen VA, Edelstein-Keshet L. Polarization and movement of keratocytes: a multiscale modelling approach. Bull Math Biol. (2006) 68:1169–211. doi: 10.1007/s11538-006-9131-7
42. Marée AF, Grieneisen VA, Hogeweg P. The Cellular Potts Model and biophysical properties of cells, tissues and morphogenesis. In: Anderson A, Chaplain M, Rejniak K. editors, Single-Cell-Based Models in Biology and Medicine, Basel: Birkhäuser. (2007), p. 107–136.
45. Käfer J, Hogeweg P, Marée AFM. Moving forward moving backward: directional sorting of chemotactic cells due to size and adhesion differences. PLoS Comput Biol. (2006) 2:e0056. doi: 10.1371/journal.pcbi.0020056
46. Segerer FJ, Thüroff F, Alberola AP, Frey E, Rädler JO. Emergence and persistence of collective cell migration on small circular micropatterns. Phys Rev Lett. (2015) 114:228102. doi: 10.1103/PhysRevLett.114.228102
47. Maiuri P, Rupprecht JF, Wieser S, Ruprecht V, Bénichou O, Carpi N, et al. Actin flows mediate a universal coupling between cell speed and cell persistence. Cell (2015) 161:374–86. doi: 10.1016/j.cell.2015.01.056
48. Garcia S, Hannezo E, Elgeti J, Joanny JF, Silberzan P, Gov NS. Physics of active jamming during collective cellular motion in a monolayer. Proc Natl Acad Sci USA. (2015) 112:15314–9. doi: 10.1073/pnas.1510973112
55. Szabó A, Varga K, Garay T, Hegedüs B, Czirók A. Invasion from a cell aggregate - The roles of active cell motion and mechanical equilibrium. Phys Biol. (2012) 9:016010. doi: 10.1088/1478-3975/9/1/016010
Keywords: cell motility, cellular Potts model, cell-cell interactions, random walk, cell adhesion
Citation: Guisoni N, Mazzitello KI and Diambra L (2018) Modeling Active Cell Movement With the Potts Model. Front. Phys. 6:61. doi: 10.3389/fphy.2018.00061
Received: 26 February 2018; Accepted: 31 May 2018;
Published: 20 June 2018.
Edited by:Francisco Monroy, Complutense University of Madrid, Spain
Reviewed by:Nir Shachna Gov, Weizmann Institute of Science, Israel
Roeland M. H. Merks, Centrum Wiskunde & Informatica, Netherlands
Copyright © 2018 Guisoni, Mazzitello and Diambra. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.