A Cellular Potts Model of the interplay of synchronization and aggregation

The model

The model uses the framework of the Cellular Potts Model (CPM; see Graner & Glazier (1992)) based on a rectangular lattice. Each cell is modeled as a group of lattice sites and is associated with an internal oscillator (“clock”). Adhesion between cells is incorporated into the model via the Hamiltonian, an energy function. Time evolution is modeled via an energy minimization procedure in space and an updating rule for the clocks. (See Fig. 1 for a sketch of the schematics of the CPM.)

More specifically, lattice sites are denoted by bold variables i = (i₁, i₂) or j = (j₁, j₂). At every time step, each lattice site i belongs to one of N cells or to the extracellular matrix (ECM). Cells are numbered with an index that ranges from 1 to N. We use the notation σ^t(i) = 0 if the site i belongs to the medium at time step t, and σ^t(i) = s if i belongs to the sth cell (s = 1, …, N). Additionally, each cell s = 1, …, N has an internal clock. This clock is a time-dependent scalar ${\theta }_s^t$. For notational simplicity, we will suppress the superscript t unless necessary.

The governing Hamiltonian is (1)${\displaystyle \mathcal{H}=\sum _{\begin{array}{c}\text{neighboring}\hfill \\ \text{lattice sites}\mathbf{i},\mathbf{j}\hfill \end{array}}\left(1-{\delta }_{\sigma \left(\mathbf{i}\right),\sigma \left(\mathbf{j}\right)}\right)f\left(\sigma \left(\mathbf{i}\right),\sigma \left(\mathbf{j}\right)\right)+\lambda \sum _{\text{cell}s}{\left(\mathrm{Area}\left(s\right)-A_{\mathrm{target}}\right)}^2.}$

The Hamiltonian encodes the compressibility of the cells and cell–cell adhesion. The second term of the Hamiltonian is a cell area constraint term. Here A(s) is the area of the cell with index s, i.e., the number of lattice sites it occupies. A_target is the target area, a fixed reference area, and λ is a parameter that encodes the compressibility of the cell: the larger λ is, the less fluctuations in size there are. (Our choices of A_target and λ are based on the work of Zhang et al. (2011); see Table 1.)

Crucially, the first term of the Hamiltonian describes cell–cell adhesion between adjacent cells. The symbol δ_ij is the Kronecker delta: δ_i,j = 1 if i = j, and δ_i,j = 0 if i ≠ j. The cell–cell adhesion term f(σ₁, σ₂) depends on the clock values of the adjacent lattice sites with indices σ₁ and σ₂ and is given by (2)${\displaystyle f\left({\sigma }_1,{\sigma }_2\right)=\left\{\begin{array}{cc}{J}_0\left(1-Jcos\left({\theta }_{{\sigma }_1}-{\theta }_{{\sigma }_2}\right)\right)\hfill & \text{if}{\sigma }_1\ne 0\text{and}{\sigma }_2\ne 0\hfill \\ J_0\hfill & \text{if}{\sigma }_1=0\text{or}{\sigma }_2=0\hfill \end{array}\right.}$

Recall that the index σ = 0 represents extracellular matrix, which has a contact energy of J₀. The parameter J encodes the effect of clock phases on adhesion. For J > 0, the function f(σ₁, σ₂) is minimized when θ_σ1 = θ_σ2, so cells of like clock phases adhere to each other the strongest (“like attracts like”). For J < 0 in contrast, f(σ₁, σ₂) is minimized when θ_σ1 = θ_σ2 + π, and so then cells of opposite clock phases adhere to each other strongest (“opposites attract”).

Time evolution is modeled via so-called spin flips. For this, a lattice site is selected randomly and it is proposed to change its index to that of a neighboring lattice site. Such a spin flip is accepted with a probability that depends on the change in the Hamiltonian function $\Delta \mathcal{H}$ it would entail. Specifically, it is given by (3)${\displaystyle \text{Prob(spin flip)}=\left\{\begin{array}{c}1\text{if}\Delta \mathcal{H}\le 0\hfill \\ e^{-\Delta \mathcal{H}/T}\text{if}\Delta \mathcal{H}>0\hfill \end{array}\right.}$

where T is the so-called temperature of the system. Higher temperatures make it more likely for spin flips to occur but also make it more likely for cells to fragment, dissolving into each other. A Monte Carlo Step (MCS) consists of the number of spin flips corresponding to the total number of sites in the lattice. Simulation time is commonly measured in MCS.

Each cell s = 1, …, N in the model has an internal oscillating clock. For each cell s = 1, …, N, we update its cell clock ${\theta }_s^t$ each MCS via (4)${\displaystyle {\theta }_s^{t+1}={\theta }_s^t+\omega \cdot \left(1+K\cdot \frac{1}{\text{\# neighbors of}s}\sum _{\text{neighbor}u\text{of}s}sin\left({\theta }_u-{\theta }_s\right)\right)}$

Note that a cell’s internal clock advances at the uniform clock speed ω, but is additionally influenced by the clocks of its neighbors in a Kuramoto-type way (Kuramoto, 1984). Here K is the clock coupling strength between neighboring cells. It controls how neighboring cells’ clocks influence each other. For K > 0, neighboring cells seek to synchronize their clock phases. For K < 0, neighboring cells seek to anti-synchronize their clocks.(The meanings of the crucial parameters J and K are summarized in Table 2.) As in other models (OKeeffe, Hong & Strogatz, 2017; Glimm & Gruszka, 2024), we note that all equations only depend on differences of clock phases, so that they are invariant under shifts in clock phase space in the form ${\theta }_s^t\to {\theta }_s^t-\omega \cdot t$. Effectively, this means we can assume that the uniform clock speed is zero and that only the terms with the K −factor in the update Eq. (4) enter into changes of the clocks at each MCS. In this sense, the choice of ω actually does not enter into the calculation apart from scaling K.

Parameters

The aim of our model is a high-level generic investigation of the interplay of adhesion and synchronization with minimal assumptions. Since it is not modeled on a concrete experimental setup, we chose to adopt the approach of Zhang et al. (2011), who used CPM simulations to investigate and validate Steinberg’s differential adhesion hypthesis (Steinberg & Takeichi, 1994). The parameter values Zhang et al. (2011) chose are based on experimental studies by Armstrong (1971) and Steinberg & Takeichi (1994). Armstrong used retinal cells from chicken embryos, Steinberg used mouse L-cells. The paraments of Zhang et al. (2011) were not the result of direct measurements, but calibrated via matching quantities such as cell size and velocity distributions of individual cells to retinal cell data (Mombach et al., 1995; Mombach & Glazier, 1996) and ensuring that cells do not fragment. Our parameter values are indicated in Table 1. The lattice length scale is approximately 2 µm per pixel, the time scale is about 10,000 MCS per hour (Zhang et al., 2011). Note that our confluency is 70%, a typical value for in vitro experiments. Our oscillator velocity ω = 0.001(clock change per MCS) corresponds to a period of about 7,000 MCS and thus very roughly on the order of one hour. This is the same order of magnitude as oscillations of Hes1 expression in somitogenesis (Kageyama, Ohtsuka & Kobayashi, 2007). We point out that in our model, only the differences between clock phases of neighboring cells matters for dynamical updates. Since all oscillators are assumed to have the same clock velocity ω, this means that simulations are actually independent of the value of ω and only depend on the effective value of the clock coordination parameter K. Accordingly, in our graphs, the color of the cells correponds to the phase shift relative to a clock that advances at a steady speed ω. We acknowledge that this independence of our simulation results on ω breaks down if we drop the assumption that every cell has the exact same speed in favor of a more realistic distribution of speeds, but this is outside the scope of the current work.

Order parameters for characterizing phases

In the Results section, we will show that our model displays qualitative different types of behaviors for different choices of the parameters J and K. This is analogous to phase states in statistical physics. Phase transitions are often described quantitatively by order parameters, approppriately defined quantities that characterize the state of the system; see e.g., Nolting (2018). We define here three order parameters for our model. In the Results section, we will use these to characterize the phases.

The first order parameter is a classical parameter quantifying synchronization due to Kuramoto (1984). It is given by (5)${\displaystyle r_{\text{global}}=\left|\frac{1}{N}\sum _{j=1}^N{e}^{i{\theta }_j}\right|}$

where N is the total number of cells in the simulation. The resulting r_global then satisfies 0 ≤ r_global ≤ 1. The parameter r_global quantifies how synchronized all cells’ clocks are at a specified MCS. If r_global = 1, all the cells’ clocks are synchronized. If r_global = 0, all the cells’ clocks are entirely unsychronized or anti-synchronized, meaning that all of their values are spread around the time clock evenly such that all values cancel each other out or there are synchronized groups of opposite phases that cancel each other out.

The second order parameter is a local version of r_global. For each cell, we compute a modified r −value via an average over its nearest neighbors (order 2). These values are then averaged over all cells: (6)${\displaystyle r_{\text{local}}=\frac{1}{N}\sum _{\text{cell}k}\left|\frac{1}{s_k}\sum _{\text{neighbor}j\text{of}k}{e}^{i{\theta }_j}\right|}$

where the first sum is taken over all cells in the simulation, s_k is the number of neighbors (order 2) of cell k and the second sum is taken over all neighbors of cell k. Note that r_global = 1 (complete global synchronization) implies r_local = 1, but not vice versa. Also note that complete randomness of phases does not typically lead to r_local = 0 because the second sum in the Eq. (6) is taken over a relatively small number, leading to large random variations. Numerically, we found that incoherence corresponds to a value of r_local of roughly 0.2; see also the discussion in the subsection “Phase Diagram”.

The third order parameter (“checkerboard parameter”) quantifies the extent to which neighboring cells’ clock phases are in opposite phase (phase difference π). We need some tolerance since neighboring cells will not be completely in opposite phase all the time. We chose a tolerance of 6.25% from total anti-synchronization between neighboring cells. More explicitly, ${\displaystyle \psi =\text{total number of pairs of neighboring lattice sites}\left(k,j\right)\text{such that}{\theta }_k-{\theta }_j\in \left[\pi -\frac{\pi }{8},\pi +\frac{\pi }{8}\right]\text{(modulo}2\pi \text{)}}$

Computational implementation

To implement the Cellular Potts Model, we used the open-source software CompuCell3D, version 4.2.5 (Swat et al., 2012) on a standard PCs (Intel(R) Core(TM) i5-6600 CPU @ 3.30 GHz, 16 GB RAM). Further analysis of the simulation data and generation of graphics were performed with Matlab. Figure 2 displays the results of 49 simulation runs, each of which took roughly 8–10 h individually.

Phase diagram

We concentrate on investigating the effects of the two parameters J and K encoding the interactions of clock synchronization and cell–cell adhesion. See Table 2 for the physical meaning of these parameters. The ranges −1 < K < 1 and −1 < J < 1 are meaningful for simulations¹. Starting with random initial conditions for cell positions and clock values, we ran simulations to determine the long-term behavior for different values of J and K. The results are summarized in the phase diagram in Fig. 2. We used movies of the simulations to confirm that the simulation duration of 250,000 MCS was sufficiently long to guarantee convergence to an obvious phase state; see also Section “Dynamics of Synchronization”. (Movies of all simulations are available online on the Zenodo repository at DOI 10.5281/zenodo.10681751).

The phase diagram in Fig. 2 clearly displays qualitatively different types of behaviors. To investigate these phases, we use the three order parameters defined in Order Section “Parameters for Characterizing Phases”. Their heat maps are shown in Fig. 3. Note that all except r_local display very sharp step-like transitions between small and large values. This allows to quantify the different phases. We identify four different phases: The first is characterized by large r_global. We call it “global synchronization.” Large values of the checkerboard parameter ψ characterize another phase, which we call “anti-synchronization.” (Specifically, we can consider ψ > 1000.) There are two more phases, which we call “incoherence” and “local synchronization”. The transition between incoherence and local synchronization is more gradual as clusters of synchronized cells become smaller with decreasing negative K. Nevertheless, we can utilize the parameter r_local to distinguish between definite local synchronization (r_local ≈ 1) and definite incoherence (r_local ≈ 0.2).

Figure 4 summarizes the phases. Global synchronization in the right half K > 0² and anti-synchronization (checkerboard pattern) in the quadrant J < 0, K ≤ 0 are both readily visibly identifiable. This is straightforward to understand: For K > 0, neighboring cells seek to synchronize; and since all cells adhere to each other, this eventually leads to global synchronization for all values of J. (The dynamics of synchronization differ by whether J is positive or negative though; we investigate this in “Dynamics of Synchronization”.) The anti-synchronized phase state is characterized by cells of opposite phases attracting each other (J < 0). Furthermore, cells seek to anti-synchronize with their neighbors (for K < 0). The resulting behavior is a checkerboard-style distribution of cell phases where cells minimize their energy by surrounding themselves with cells of the opposite phase.

The other two phases (local synchronization (synchronized spatial clusters) and incoherence (quasi-random spatial distribution of phases)) both occur in the quadrant J ≥ 0, K ≤ 0. For very negative K and small J, we have incoherence. In the case large J and K close to 0, cells tend to sort by phases and clusters de-synchronize sufficiently slowly that persistent locally synchronized clusters form.

Movies show that global synchronization and anti-synchronization are essentially static distributions, where the cells’ positions and shapes fluctuate, but the overall spatial distribution of clock phases stays the same. Incoherence and local synchronization are dynamic phases, where cells or cell clusters in different clock phases constantly move relative to each other, yielding a behavior in which each snapshot in time is different, but characteristic of the typical distribution.

To further characterize the incoherence phase, we compared the resulting distributions of clock values to the distribution obtained by chance. The results are summarized in Fig. 5. Note that the distributions of clock values for the incoherent phase is essentially indistinguishable to a uniform random distribution on the interval [0, 2π). There is a marked difference to the anti-synchronization (checkerboard) phase, with two peaks at the distance π. Interestingly, the “local synchronization” case also gives a distribution indistinguishable from chance. This is because in the representation of Fig. 5, all spatial information is lost, so spatially completely mixed clock phases and spatially clustered clock phases give similar distributions.

Dynamics of synchronization

For positive values of K, cells’ clocks eventually globally synchronize as shown in Fig. 2. When we investigated the dynamical paths to synchrony however, we found substantial qualitative differences depending on the sign of the parameter J (positive, negative or zero). Figure 6 illustrates this with K = 1 and three different values of J. To quantify the initial progress of synchronization for different values of J, Fig. 7A shows the order parameter r_global over time (MCS) up to 20,000 MCS. There is a clear hierarchy of J-values –the speed of synchronization is highest for the most negative value of J and decreases with increasing J with J =  + 0.95 corresponding to the slowest rate of progress. This can be seen also by the first two columns in Fig. 6. Why is this? Again, Fig. 6 provides an important insight: For negative J, “opposites attract” and hence at t = 10, 000 MCS, it is clearly visible that for J =  − 0.95, there are large regions of synchronized cells with cells of the opposite phase interspersed. (For example regions of orange cells with interspersed small blue cells, some slightly fragmented.) This proximity of cells of opposite clock phases early on leads to rapid synchronization. In contrast, this is not the case for J =  + 0.6333 (positive J, so “like attracts like”) or J = 0 (clock phase has no influence on adhesion). One also observes that synchronized regions have more gradual transitions for J = 0 than for J =  + 0.6333. For instance for J =  + 0.6333, red and green cells tend to be separated by much smaller buffers of yellow and orange cells than is the case for J = 0.

These observations can be quantified via the local parameter r_local, shown in Fig. 7C: Here, the order is reversed and initially, the smaller J, the slower the rate of local synchronization. For negative J (“opposites attract”), the local synchronization even decays before increasing again.

The long-term rate of synchronization is plotted in Fig. 7B. The relationship between order of J −values and speed of synchronization seen in Fig. 7A is not preserved in the long run. Some of this may be due to the small sample size (n = 10). However, it is very clear that synchronization progresses most slowly for J = 0, the case where clock phase does not influence adhesion. At 250,000 MCS, complete synchronization had not been reached. (Illustrated also in Fig. 6 and Fig. 4.) This is in contrast to positive or negative values of J, where complete synchronization was achieved at that time. (Even for J = 0, synchronization was eventually achieved if the simulation was allowed to keep running; see time t = 400, 000 MCS in Fig. 6.)

Why is synchronization especially slow for J = 0? An intuitive explanation is as follows: Synchronization proceeds especially fast (for positive K) if the clock phase distribution has sharp spatial gradients, i.e., if many cells of very different clock phases are neighbors. In contrast, configurations with shallow clock gradients (gradual transitions of clock values) exhibit slower synchronization. With this principle, it is clear that for negative J, “opposites attract” and one indeed gets fields of cells with sharp clock gradients which synchronize quickly, as we noted before. But also the case of positive J (“like attracts like”) creates sharp gradients, but with a different mechanism: Early on, synchronized clusters form as cells with similar phases effectively move towards each other. Because of cell–cell adhesion, these clusters have relatively sharply defined edges (Fig. 6), meaning sharp clock gradients, again leading to faster synchronization. In contrast, for J = 0, clock gradients are much more gradual, allowing regions of different clock phases to persist longer.