Theory

General Cellular Potts Model (CPM)

The Cellular Potts Model (CPM) is a lattice-based computational framework used to simulate the collective behavior of cells, such as their movement, growth, and interactions. It is an extension of the classical Potts model from statistical mechanics and is particularly effective in modeling biological processes where cell shape and adhesion are critical factors.

Lattice and Cell Representation

In the CPM, the domain is represented as a discrete lattice (typically a 2D grid) where each site \((i, j)\) is assigned a cell identifier \(\sigma(i, j)\). The value of \(\sigma(i, j)\) can be:

\[ \sigma(i, j) = \begin{cases} 0 & \text{if the site is part of the medium (empty space)},\\ k & \text{if the site belongs to cell } k \text{ (where } k \geq 1\text{)}. \end{cases} \]

Cells are therefore clusters of lattice sites with the same \(\sigma\) value, and the medium is denoted by \(\sigma = 0\).

Hamiltonian and Energy Terms

The behavior of cells in the CPM is governed by an energy function (Hamiltonian). The system evolves by attempting to minimize this Hamiltonian through stochastic processes. The Hamiltonian typically includes several energy terms to capture different physical properties of cells.

Interaction Energy

The interaction energy models cell adhesion and is defined as:

\[ H_{\text{interaction}} = \sum_{i} \sum_{j} \sum_{\sigma' \in \text{neighbors}} J_{\sigma(i,j), \sigma'} (1 - \delta_{\sigma(i,j), \sigma'}), \]

where:

• \(J_{\sigma(i,j), \sigma'}\) is the interaction strength between cells \(\sigma(i,j)\) and \(\sigma'\).
• \(\delta_{\sigma(i,j), \sigma'}\) is the Kronecker delta, which equals 1 if \(\sigma(i,j) = \sigma'\) and 0 otherwise.

This term penalizes boundaries between cells, promoting adhesion between similar cell types while resisting adhesion between different types.

Surface Energy

To maintain cell shape, a surface energy constraint is imposed:

\[ H_{\text{surface}} = \sum_{i}^N \lambda_{\sigma(i)}^{S} (a_i - A_{\sigma(i)})^2, \]

where:

• \(a_i\) is the actual surface area of cell \(i\).
• \(A_{\sigma(i)}\) is the target surface area for the cell type \(\sigma(i)\).
• \(\lambda_{\sigma(i)}^{S}\) is a parameter controlling the strength of the surface constraint.

This term ensures that cells maintain a size close to their target surface area.

Monte Carlo Dynamics

The CPM evolves by performing Monte Carlo simulations to minimize the Hamiltonian. In each Monte Carlo step, a site \((i, j)\) on the lattice is randomly chosen, and a copy attempt is made to change its identifier \(\sigma(i, j)\) to that of a neighboring site. The change in energy \(\Delta H\) due to the proposed copy is calculated. The copy is accepted with a probability:

\[ p = \begin{cases} 1 & \text{if } \Delta H < 0,\\ e^{-\Delta H / kT} & \text{otherwise}, \end{cases} \]

where:

• \(k\) is the Boltzmann constant.
• \(T\) is the system temperature, controlling the level of randomness.

This mechanism allows the system to explore configurations that lower the overall energy while occasionally accepting higher-energy changes, facilitating the simulation of dynamic cellular processes.

Time Unit: Monte Carlo Sweep (MCS)

A single Monte Carlo Sweep (MCS) consists of as many copy attempts as there are lattice sites. This serves as the basic time unit in CPM simulations, representing one iteration of the system's evolution.

By iterating over many MCS, the cells gradually change their configuration, adhering to physical constraints and interaction rules defined by the Hamiltonian.

My Potts Model

For \( N \) cells of type \( \sigma(i) \), the Hamiltonian is defined as follows:

\[ H = H_{\text{interaction}} + H_{\text{surface}} + H_{\text{length}} + H_{\text{chemotaxis}} + \epsilon_{\text{connectivity}} \]

with:

\[ H_{\text{interaction}} = \sum_{i} \sum_{j} \sum_{\sigma'\in\text{ neighbors}} J_{\sigma(i,j), \sigma'} (1 - \delta_{\sigma(i,j),\sigma'}) \]

\[ H_{\text{surface}} = \sum_{i}^N \left( \lambda_{\sigma(i)}^{S} (a_i - A_{\sigma(i)})^2 \right) \]

\[ H_{\text{length}} = \sum_{i}^N \left( \lambda_{\sigma(i)}^{L} (l_i - L_{\sigma(i)})^2 \right) \]

\[ \Delta H_{\text{chemotaxis}} = \sum_{i}\sum_{j} \left( \lambda_{\sigma(i,j)}^{C} \Delta C(i,j) \right) \]

with \( J \) the interaction matrix, \( \lambda_{\sigma(i)}^{S} \) and \( \lambda_{\sigma(i)}^{L} \) the surface and length parameters, \( a_i \) and \( l_i \) the surface and length of the cell \( i \), \( A_{\sigma(i)} \) and \( L_{\sigma(i)} \) the target surface and length of the cell type \( \sigma(i) \), \( \lambda_{\sigma(i)}^{C} \) the chemotaxis parameter, and \( \nabla C(i,j) \) the gradient of the chemoattractant at position \( (i,j) \).
The connectivity term \( \epsilon_{\text{connectivity}} \) is a penalty term that ensures that the cells do not break apart ( \( \epsilon_{\text{connectivity}} \geq 20 kT \)).

The Chemotaxis term is not calculated for both sites but just added to the general \( \Delta H \). It is using another lattice with the same size as the lattice of the cells. This lattice lets us have a concentration of chemoattractant at each position. The evolution of chemoattractant is calculated with the following equation:

\[ \frac{\partial C(x,y,t)}{\partial t} = D \nabla^2 C(x,y,t) + \eta(x,y) (1 - C(x,y,t)) - \gamma(x,y) C(x,y,t) \]

with \( D \) the diffusion coefficient, \( \eta(x,y) \) the secretion rate at location \( (x,y) \), \( \gamma \) the degradation rate at location \( (x,y) \), and \( \nabla^2 C(x,y,t) \) the Laplacian of the concentration at location \( (x,y) \).

This will be solved with a finite difference method:

\[ C(i,j,t+1) = C(i,j,t) + \Delta t \left[ D \nabla^2 C(i,j,t) + \eta(1-\delta(\sigma(i,j),0)) - \gamma\delta(\sigma(i,j),0) C(i,j,t) \right] \]

with

\[ \nabla^2 C(i,j,t) = C(i+1,j,t) + C(i-1,j,t) + C(i,j+1,t) + C(i,j-1,t) - 4 C(i,j,t) \]

and with the secretion rate \( \eta \) and the degradation rate \( \gamma \) defined as follows:

\[ \begin{cases} \eta(x,y) = \eta_0 \cdot \delta(\sigma(i,j),1) \\ \gamma(x,y) = \gamma_0 \cdot \delta(\sigma(i,j),0) \end{cases} \]

The Metropolis algorithm is used to minimize the energy of the system by moving the cells on the lattice. At each Monte Carlo step, a random cell is selected, and a random move is proposed. The move is accepted if the change in energy \( \Delta H \) is negative or with a probability \( p = e^{-\Delta H / kT} \):

\[ p = \begin{cases} 1 & \text{if } \Delta H < 0 \\ e^{-\Delta H / kT} & \text{otherwise} \end{cases} \]

where \( \Delta H \) is the change in energy due to the move, \( k \) is the Boltzmann constant, and \( T \) is the temperature of the system.

Time unit

The time unit is a Monte Carlo Sweep (MCS), which corresponds with as many copy attempts as there are lattice sites.

Inertia Tensor

The inertia tensor will be calculated for each cell. This will allow us to have the orientation of the cells and especially the elongation of the cell. This will be used for a new term in the Hamiltonian.

The Inertia tensor is defined as:

\[ I = \begin{pmatrix} I_{xx} & - I_{xy} \\ - I_{yx} & I_{yy} \end{pmatrix} \]

with the diagonal terms:

\( I_{xx} = \sum_{i,j} m_{ij} (y_i - y_c)^2 \) and \( I_{yy} = \sum_{i,j} m_{ij}(x_j - x_c)^2 \)

and the off-diagonal terms (without the minus sign):

\( I_{xy} = I_{yx} = \sum_{i,j} m_{ij} (x_i - x_c)(y_j - y_c) \)

Connectivity Constraint

To check whether an index copy into site \( x \) will change local connectivity, we count how many of site \( x \)'s neighbors \( x_i \) have the same index as \( x \) while the next location clockwise around \( x \) has an index different from that of \( x \). If this quantity:

\[ \sum_i \delta_{\sigma(x),\sigma(x_i)}(2 - \delta_{\sigma(x),\sigma(x_{i+1})} - \delta_{\sigma(x),\sigma(x_{i-1})}) > 2 \]

changing the site will destroy the local connectivity.