Potts and Lattice
Potts Algorithm: How It Works
- Learning Objectives:
Understand how the Potts algorithm determines cell shape by controlling their pixels
The first section of the .xml file defines the global parameters of the lattice and the simulation. Your file may look much simpler than this.
<Potts>
<Dimensions x="101" y="101" z="1"/>
<Anneal>0</Anneal>
<Steps>1000</Steps>
<FluctuationAmplitude>5</FluctuationAmplitude>
<Flip2DimRatio>1</Flip2DimRatio>
<Boundary_y>Periodic</Boundary_y>
<Boundary_x>Periodic</Boundary_x>
<NeighborOrder>2</NeighborOrder>
<DebugOutputFrequency>20</DebugOutputFrequency>
<RandomSeed>167473</RandomSeed>
<EnergyFunctionCalculator Type="Statistics">
<OutputFileName Frequency="10">statData.txt</OutputFileName>
<OutputCoreFileNameSpinFlips Frequency="1" GatherResults="" OutputAccepted="" OutputRejected="" OutputTotal=""/>
</EnergyFunctionCalculator>
</Potts>
This section appears at the beginning of the configuration file.
The line reading <Dimensions x="101" y="101" z="1"/> declares the dimensions of the
lattice to be 101 by 101 by 1 pixels.
Since z is small, this lattice is two-dimensional.
Lattice sites are 0-indexed, which means that we count them from 0 to 100 when
accessing them in Python.
<Steps>1000</Steps> tells CompuCell how
long the simulation lasts in Monte Carlos Steps (MCS). After executing this number of steps, CompuCell can simulate with zero temperature for an additional
period. In our case, it will run for <Anneal>10</Anneal> extra steps.
Fluctuation Amplitude (Temperature)
The FluctuationAmplitude/Temperature parameter determines the intrinsic fluctuation or
motility of each cell membrane. (Either key word works in XML)
Try setting a high Temperature value in XML.
You should see much more activity from pixel changes.
Conversely, a low Temperature will decrease cell membrane activity.
Note
FluctuationAmplitude is a Temperature
parameter in classical GGH Model formulation. We have decided to use
FluctuationAmplitude term instead of temperature because using the word
Temperature to describe the intrinsic motility of cell membranes was quite
confusing.
In the above example, fluctuation amplitude applies to all cells in the simulation. To define fluctuation amplitude separately for each cell type, we use the following syntax. Each value must be non-negative.
<FluctuationAmplitude>
<FluctuationAmplitudeParameters CellType="Condensing" FluctuationAmplitude="10"/>
<FluctuationAmplitudeParameters CellType="NonCondensing" FluctuationAmplitude="5"/>
</FluctuationAmplitude>
When CompuCell3D encounters the expanded definition of FluctuationAmplitude,
it will use it in place of a global definition like this one:
<FluctuationAmplitude>5</FluctuationAmplitude>
Alternatively, you can use Python to set the fluctuation amplitude individually for each cell:
for cell in self.cellList:
if cell.type==1:
cell.fluctAmpl=20
When determining which value of fluctuation amplitude to use, CompuCell
prioritizes Python definitions. Otherwise, if fluctAmpl was not set by Python, it will try to use the CC3DML for fluctuation amplitude by cell types.
Lastly, it will resort to a globally defined fluctuation amplitude
(Temperature). If none of these are defined, the default value is used. Thus, it is perfectly fine to use a combination of these techniques.
Note
Default Value of Fluctuation Amplitude
For XML-based simulations, the default value of FluctuationAmplitude/Temperature is 0. For Python-only (PyCoreSpecs) simulations, the default value of fluctuation_amplitude is 10.
Note that a value of 0 does not turn off activity; there will still be slight fluctuations in cell membranes.
In the Glazier-Graner-Hogeweg (GGH) Model, the fluctuation amplitude is determined by taking into account the fluctuation amplitude of a “source” (expanding) cell and a “destination” cell (the one that will be overwritten).
Currently, CompuCell3D supports functions used to calculate resultant fluctuation amplitude (those functions take as argument fluctuation amplitude of “source” and “destination” cells and return fluctuation amplitude that is used in
calculation of pixel-copy acceptance). The 3 functions are Min, Max, and
ArithmeticAverage and we can set them using the following option of the
Potts section:
<Potts>
<FluctuationAmplitudeFunctionName>Min</FluctuationAmplitudeFunctionName>
…
</Potts>
By default, we use the Min function. Notice that if you use the global
fluctuation amplitude definition Temperature, it does not really matter
which function you use. The differences arise when “source” and
“destination” cells have different fluctuation amplitudes.
The above concepts are best illustrated by the following example:
<Potts>
<Dimensions x="100" y="100" z="1"/>
<Steps>10000</Steps>
<FluctuationAmplitude>5</FluctuationAmplitude>
<FluctuationAmplitudeFunctionName>ArithmeticAverage</FluctuationAmplitudeFunctionName>
<NeighborOrder>2</NeighborOrder>
</Potts>
Where in the CC3DML section we define global fluctuation amplitude and
we also use ArithmeticAverage function to determine resultant
fluctuation amplitude for the pixel copy.
Try this Python script to see how fluctuation amplitude affects
the membranes of cells.
This code assigns a different fluctAmpl value depending on
which of 4 quadrants each cell is located in.
class FluctuationAmplitude(SteppableBasePy):
def __init__(self, _simulator, _frequency=1):
SteppableBasePy.__init__(self, _simulator, _frequency)
self.quarters = [[0, 0, 50, 50], [0, 50, 50, 100], [50, 50, 100, 100], [50, 0, 100, 50]]
self.steppableCallCounter = 0
def step(self, mcs):
quarterIndex = self.steppableCallCounter % 4
quarter = self.quarters[quarterIndex]
for cell in self.cellList:
if cell.xCOM >= quarter[0] and cell.yCOM >= quarter[1] and cell.xCOM < quarter[2] and cell.yCOM < quarter[3]:
cell.fluctAmpl = 50
else:
# this means CompuCell3D will use globally defined FluctuationAmplitude
cell.fluctAmpl = -1
self.steppableCallCounter += 1
Similarly, fluctuation_amplitude can be set in a Python-only simulation:
from cc3d import CompuCellSetup
from cc3d.core.PyCoreSpecs import Metadata, PottsCore
spec_potts = PottsCore()
spec_potts.dim_x, spec_potts.dim_y = 100, 100
spec_potts.steps = 100000
spec_potts.neighbor_order = 2
spec_potts.fluctuation_amplitude = 50
OR
from cc3d import CompuCellSetup
from cc3d.core.PyCoreSpecs import Metadata, PottsCore
CompuCellSetup.register_specs(PottsCore(dim_x=dim_x,
dim_y=dim_y,
steps=100000,
neighbor_order=2,
boundary_x="Periodic",
boundary_y="Periodic",
fluctuation_amplitude=50))
Remember, negative values of fluctuationAmplitude are ignored.
Here, cell.fluctAmpl = -1 is a hint to CC3D to use fluctuation amplitude defined in the
CC3DML.
Let us revisit our original example of the Potts section CC3DML:
<Potts>
<Dimensions x="101" y="101" z="1"/>
<Anneal>0</Anneal>
<Steps>1000</Steps>
<FluctuationAmplitude>5</FluctuationAmplitude>
<Flip2DimRatio>1</Flip2DimRatio>
<Boundary_y>Periodic</Boundary_y>
<Boundary_x>Periodic</Boundary_x>
<NeighborOrder>2</NeighborOrder>
<DebugOutputFrequency>20</DebugOutputFrequency>
<RandomSeed>167473</RandomSeed>
<EnergyFunctionCalculator Type="Statistics">
<OutputFileName Frequency="10">statData.txt</OutputFileName>
<OutputCoreFileNameSpinFlips Frequency="1" GatherResults="" OutputAccepted="" OutputRejected="" OutputTotal=""/>
</EnergyFunctionCalculator>
</Potts>
Based on our discussion about the difference between pixel-flip attempts and
MCS (see “Introduction to CompuCell3D”), we can specify how many pixel
copies should be attempted in every MCS. We specify this number
indirectly by specifying the Flip2DimRatio by using
<Flip2DimRatio>1</Flip2DimRatio>
which tells CompuCell that it should
make 1 times number of lattice sites attempts per MCS – in our case one MCS
is 101x101x1 pixel-copy attempts. To set 2.5 x 101 x 101 x 1 pixel-copy
attempts per MCS you would write:
<Flip2DimRatio>2.5</Flip2DimRatio>
The line beginning with <NeighborOrder>2</NeighborOrder> specifies the neighbor order.
Neighbor order controls how many nearby pixels the Potts algorithm will check
each time it needs to do a pixel copy attempt.
Think of the neighbors as a circular area around each pixel.
If you set a higher neighbor order, you may have smoother cells but less performance.
In the previous example, the pixel neighbors are ranked according to their distance from a reference pixel (i.e. the one
you are measuring a distance from). Thus, we can group the 1st, 2nd, and 3rd nearest neighbors for every pixel in the lattice. Using 1st nearest neighbor
interactions may cause unwanted artifacts due to lattice anisotropy. The longer the interaction range,
(i.e. 2nd, 3rd or higher NeighborOrder), the more isotropic the
simulation and the slower it runs. In addition, if the interaction range
is comparable to the cell size, you may generate unexpected effects,
since non-adjacent cells will contact each other.
On a hex lattice, those problems seem to be less severe and there 1st or 2nd nearest neighbor usually are sufficient.
Ranking of pixel neighbors on square 2D lattice
The Potts section also contains tags called <Boundary_y> and
<Boundary_x>. These tags impose boundary conditions on the lattice. In
this case, the x and y axes are periodic.
Periodic Boundary Conditions: cause the edges of the simulation area to “wrap around.” For example, a pixel at (x=0 , y=1, z=1)
will neighbor the pixel at (x=100, y=1, z=1). We recommend using periodic boundaries when you want to simulate a large area of tissue while keeping your lattice small.
<Boundary_x>Periodic</Boundary_x>
‘NoFlux’ Boundary Conditions: is the opposite of periodic, so the lattice remains a finite area. This is the default.
<Boundary_x>NoFlux</Boundary_x>
Boundary conditions are independent in each XYZ direction, so you can specify any combination of them you like.
See Boundary Conditions for more details.
DebugOutputFrequency: is used to tell CompuCell3D how often it should output text information about the status of the simulation. This tag is optional.
RandomSeed: is used to initialize the random number generator. You do not need this tag unless you want every simulation to behave exactly the same, which is not recommended. See Stochasticity and RandomSeed for more details.
EnergyFunctionCalculator: allows you to output statistical data, such as the changes in energy from the simulation, to text files for further analysis. See How to Output Energy Changes for more details.
One option of the Potts section that we have not used here is the ability to customize acceptance function for Metropolis algorithm:
<Offset>-0.1</Offset>
<KBoltzman>1.2</KBoltzman>
This ensures that pixel copy attempts that increase the energy of the system are accepted with probability
where \(δ\) and \(k\) are specified by Offset and KBoltzman tags, respectively.
By default, \(δ=0\) and \(k=1\). (That is, Offset is 0 and KBoltzman is 1).
As an alternative to the exponential acceptance function, you may use a simplified version, which is essentially 1 order of expansion of the exponential:
To be able to use this function, all you need to do is to add the following line in the Potts section:
<AcceptanceFunctionName>FirstOrderExpansion</AcceptanceFunctionName>
Lattice Type
Early versions of CompuCell3D allowed users to use only square lattice. Most recent versions allow the simulation to be run on hexagonal lattice as well.
Full description of hexagonal lattice including detailed derivations can be found in “Introduction to Hexagonal Lattices” available from http://www.compucell3d.org/BinDoc/cc3d_binaries/Manuals/HexagonalLattice.pdf
To enable hexagonal lattice you need to put
<LatticeType>Hexagonal</LatticeType>
in the Potts section of the CC3DML configuration file.
There are few things to be aware of when using hexagonal lattice.
In 2D your pixels are hexagons but in 3D the voxels are rhombic dodecahedrons.
It is particularly important to realize that surface or perimeter of the pixel
(depending whether in 2D or 3D) is different than in the case of square
pixel. The way CompuCell3D hex lattice implementation was done was that
the volume of the pixel was constrained to be 1 regardless of the
lattice type.
There is also one to one correspondence between pixels of the square
lattice and pixels of the hex lattice. Consequently, we can come up with
transformation equations which give positions of hex pixels as a
function of square lattice pixel position:
Based on the above facts one can work out how unit length and unit surface transform to the hex lattice. The conversion factors are given below:
For the 2D case, assuming that each pixel has unit volume, we get:
where \(S_{hex-unit}\) denotes length of the hexagon and \(L_{hex-unit}\) denotes a distance between
centers of the hexagons. Notice that unit surface in 2D is simply a
length of the hexagon side and surface area of the hexagon with side a
is:
In 3D we can derive the corresponding unit quantities starting with the formulae for volume and surface of rhombic dodecahedron (12 hedra)
where a denotes length of dodecahedron edge.
Constraining the volume to be 1 we get:
and thus unit surface is given by:
and unit length by: