Cellular Potts model of cell sorting on square, triangular, and disordered lattices
The cellular Potts model of a confluent tissue whose cells of several kinds sort by differential adhesion.
- Name:
cellular_potts - Field:
biology/physiology - Formalism: lattice dynamics
- References: Nemati and de Graaf 2024, arXiv 2404.09055; Graner and Glazier 1992
Select this system
from skcomplex.systems import create_simulation, create_system, describe
system = create_system('cellular_potts')
entry = describe('cellular_potts')
trajectory = create_simulation('cellular_potts').prepare().run().trajectory_
Factory
@register_system(
field='biology/physiology',
tags=('cellular potts', 'tissue', 'cell sorting', 'differential adhesion', 'lattice', 'monte carlo'),
references=('Nemati and de Graaf 2024, arXiv 2404.09055', 'Graner and Glazier 1992'),
defaults={'until': 2000, 'dt': 1},
observables=_cellular_potts_observables,
)
def cellular_potts(
shape: Sequence[int] = (40, 40),
*,
lattice: str = 'square',
cell_shape: Sequence[int] = (5, 8),
adhesion: npt.ArrayLike = DEFAULT_ADHESION,
temperature: float = 2.0,
area_strength: float = 1.0,
perimeter_strength: float = 0.0,
target_perimeter: float = 0.0,
hits: int = 1,
random_state: int | np.random.Generator | None = 0,
) -> System:
"""Return the cellular Potts model of a confluent tissue whose cells of several kinds sort by differential adhesion.
Every site of a periodic lattice belongs to one cell, and the cells start as rectangles of `cell_shape` sites, whose
area is the target area, with kinds in equal shares at random. Label copies sample the energy
`sum_<i,j> w_ij J(t_i, t_j) [s_i != s_j] + lambda sum_s (a_s - A)^2 + mu sum_s (p_s - P)^2` of Nemati and de Graaf
(2024), Eq. 3, extended to the adhesion `J` between kinds of Graner and Glazier (1992). The square lattice reads the
Moore neighbourhood and the triangular lattice the six nearest neighbours, both with unit weights, and the
disordered lattice is the periodic Voronoi tessellation of one jittered point per unit area, whose sites have the
area of their Voronoi cell and touch along sides of weight `w_ij`, the length of the side over the mean length. A
heterotypic adhesion above the mean of the homotypic ones raises the share of contacts between cells of one kind.
Args:
shape:
The number of rows and columns of the lattice, or of points of the disordered lattice, multiples of the
cell shape.
lattice:
`'square'`, `'triangular'`, or `'disordered'`.
cell_shape:
The rows and columns of the initial rectangle of every cell, whose area is the target area.
adhesion:
The adhesion energy per unit contact between kinds, a symmetric matrix with one row per kind.
temperature:
The temperature of the label copies.
area_strength:
The strength `lambda` of the area constraint.
perimeter_strength:
The strength `mu` of the perimeter constraint.
target_perimeter:
The target perimeter `P`, the weighted number of contacts of a cell with other cells.
hits:
The number of copies proposed per site and sweep.
random_state:
The seed or generator of the kinds of the cells and of the disordered lattice.
Returns:
system:
The system, with the states `tissue.cell`.
Raises:
ValidationError:
The lattice is unknown, or the shape is not a multiple of the cell shape.
Examples:
>>> from skcomplex.systems.biology.physiology import cellular_potts
>>> system = cellular_potts((10, 16))
>>> system.interactions[0].kinds
(0, 0, 1, 1)
>>> cellular_potts((10, 16), lattice='disordered').elements[0].space.graph.number_of_nodes()
160
"""
rows, columns = (int(size) for size in shape)
height, width = (int(size) for size in cell_shape)
if lattice not in LATTICES:
error_msg = f'The lattice of the cellular Potts model is one of {", ".join(LATTICES)}, but got {lattice!r}.'
raise ValidationError(error_msg)
if rows % height or columns % width:
error_msg = f'The shape {(rows, columns)} of the tissue is not a multiple of the cell shape {(height, width)}.'
raise ValidationError(error_msg)
rng = np.random.default_rng(random_state)
count = (rows // height) * (columns // width)
kinds = np.asarray(adhesion, dtype=float).shape[0]
kind_of = tuple(int(kind) for kind in rng.permutation(np.arange(count) % kinds))
if lattice == 'disordered':
graph = _disordered_lattice((rows, columns), rng)
x, y = np.asarray([point for _, point in graph.nodes(data='pos')]).T
space: LatticeSpace | NetworkSpace = NetworkSpace(graph)
else:
y, x = np.indices((rows, columns)).reshape(2, -1).astype(float)
space = LatticeSpace((rows, columns), geometry=lattice)
cells = (y // height).astype(int) * (columns // width) + (x // width).astype(int)
tissue = LatticeElement('tissue', element=TissueSiteElement, space=space, values={'cell': cells.astype(float)})
sampling = CellAdhesionInteraction(
'adhesion',
target='tissue',
kinds=kind_of,
adhesion=adhesion,
temperature=temperature,
area_strength=area_strength,
target_area=float(height * width),
perimeter_strength=perimeter_strength,
target_perimeter=target_perimeter,
hits=hits,
)
return System(elements=[tissue], interactions=[sampling], time=DiscreteTime())
Components
Elements:
TissueSiteElement, fromElement
Interactions:
CellAdhesionInteraction, fromSamplingInteraction
class TissueSiteElement(Element):
"""Site of a lattice tissue, which holds the label of the cell that covers it.
The sites of one label form a cell of the cellular Potts model, and `CellAdhesionInteraction` moves the boundaries
between cells by copying the label of a neighbouring site.
Args:
label:
The label of the site or the lattice group.
visual:
How the site is drawn, or `None` for `default_visual_`.
cell:
The initial label of the cell that covers the site, from 0.
Attributes:
default_visual_:
The visual of the site when it sets none.
Examples:
>>> from skcomplex.systems.biology.physiology import TissueSiteElement
>>> TissueSiteElement('site', cell=3).cell
3
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
cell: Annotated[float, StateSpec(domain=NN)] = 0.0,
) -> None:
super().__init__(label, visual=visual)
self.cell = cell
class CellAdhesionInteraction(SamplingInteraction):
"""Cellular Potts sampling of a tissue, whose cells move their boundaries under adhesion and area constraints.
The energy of the Graner-Glazier model is `sum_<i,j> w_ij J(t(s_i), t(s_j)) [s_i != s_j] + lambda sum_s (a_s -
A)^2 + mu sum_s (p_s - P)^2`, where `s_i` is the cell of site `i`, `t(s)` the kind of cell `s`, `J` the adhesion
energy per unit contact between kinds, `w_ij` the weight of the contact of two sites, `a_s` the area of a cell, and
`p_s` its perimeter, the weighted number of its contacts with other cells. Every sweep proposes `hits` copies per
site of the cell of a random neighbour into a random site, accepted by the Metropolis rule at the temperature. A
heterotypic adhesion above the mean of the homotypic ones sorts the kinds, the differential adhesion of Steinberg.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the group of sites.
kinds:
The kind of every cell, by label.
adhesion:
The adhesion energy per unit contact between cells of every pair of kinds, a symmetric matrix.
temperature:
The temperature of the copies, the amplitude of the fluctuations of the cell boundaries.
area_strength:
The strength `lambda` of the area constraint.
target_area:
The target area `A` of every cell.
perimeter_strength:
The strength `mu` of the perimeter constraint.
target_perimeter:
The target perimeter `P` of every cell.
state:
The name of the state that holds the cell of a site.
rule:
The acceptance, `'metropolis'` or `'glauber'`.
neighborhood:
The neighbours of a site, or `None` for the default neighbourhood of the space.
hits:
The number of copies proposed per site and sweep.
Raises:
ValidationError:
The adhesion is not a symmetric square matrix, or a kind is not one of its rows.
Examples:
>>> from skcomplex.systems.biology.physiology import CellAdhesionInteraction
>>> adhesion = CellAdhesionInteraction('adhesion', target='tissue', kinds=(0, 1), adhesion=((1, 2), (2, 1)))
>>> adhesion.objective.contact(np.array([0, 0]), np.array([0, 1])).tolist()
[1.0, 2.0]
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
kinds: Sequence[int],
adhesion: npt.ArrayLike,
temperature: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
area_strength: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
target_area: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 40.0,
perimeter_strength: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
target_perimeter: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
state: str = 'cell',
rule: str = 'metropolis',
neighborhood: BaseNeighborhood | None = None,
hits: int = 1,
) -> None:
super().__init__(
label,
target=target,
state=state,
objective=RegionScoreSpec(self._get_contact, self._get_area, self._get_perimeter, hits),
move='copy',
rule=rule,
neighborhood=neighborhood,
)
self.kinds = kinds
self.adhesion = adhesion
self.temperature = temperature
self.area_strength = area_strength
self.target_area = target_area
self.perimeter_strength = perimeter_strength
self.target_perimeter = target_perimeter
self.hits = hits
def _get_intensity(self, namespace: SimpleNamespace) -> object:
"""Return the intensity of the acceptance, one over the temperature."""
return 1 / SR(self._get_parameter_symbol('temperature', namespace))
def _get_contact(self, first: Indices, second: Indices) -> Array:
"""Return the adhesion energy per unit contact between the cells of two arrays of labels."""
kinds = np.asarray(self.kinds, dtype=np.intp)
return np.asarray(self.adhesion, dtype=float)[kinds[first], kinds[second]]
def _get_area(self, labels: Indices, area_strength: float, target_area: float) -> tuple[float, float]:
"""Return the strength and the target of the area constraint, the same for every cell."""
del labels
return area_strength, target_area
def _get_perimeter(
self,
labels: Indices,
perimeter_strength: float,
target_perimeter: float,
) -> tuple[float, float]:
"""Return the strength and the target of the perimeter constraint, the same for every cell."""
del labels
return perimeter_strength, target_perimeter
def _check(self) -> None:
"""Raise an error for an adhesion that is not a symmetric square matrix, or a kind that is not one of its rows.
Raises:
ValidationError:
The adhesion is not a symmetric square matrix, or a kind is not one of its rows.
"""
super()._check()
adhesion = np.asarray(self.adhesion, dtype=float)
if (
adhesion.ndim != MATRIX_AXES
or adhesion.shape[0] != adhesion.shape[1]
or not np.allclose(adhesion, adhesion.T)
):
error_msg = f'The adhesion of {self.label!r} is a symmetric matrix over the kinds, but got {adhesion!r}.'
raise ValidationError(error_msg)
kinds = np.asarray(self.kinds)
if kinds.ndim != 1 or (kinds.size and (kinds.min() < 0 or kinds.max() >= adhesion.shape[0])):
error_msg = f'The kinds of {self.label!r} are rows of the adhesion, from 0 to {adhesion.shape[0] - 1}.'
raise ValidationError(error_msg)