Skip to content

Commit

Permalink
Merge pull request #27 from CompuCell3D/textbook
Browse files Browse the repository at this point in the history
Rewrote Steppable, Potts, Contact Plugin, & Cell Sorting
  • Loading branch information
maciekswat committed Jan 6, 2024
2 parents 1de0406 + 599c4b1 commit 6c72b00
Show file tree
Hide file tree
Showing 45 changed files with 872 additions and 512 deletions.
103 changes: 30 additions & 73 deletions docs/SteppableBasePy_class.rst
@@ -1,92 +1,49 @@
SteppableBasePy class
=====================
What is a Steppable? (SteppableBasePy class)
===================================================

In the example above you may wonder how it is possible that it is
sufficient to type:
``SteppableBasePy`` has built-in **functions that are called automatically** during the simulation.
The most important functions are ``start`` and ``step``.

.. code-block:: python
for cell in self.cell_list:
to iterate over a list of all cells in the simulation. Where does
``self.cell_list`` come from and how it accesses/stores information about all
cells? The full answer to this question is beyond the scope of this
manual so we will give you only a hint what happens here. The
``self.cell_list`` is a member of a ``SteppableBasePy`` class. All CC3D Python
steppable inherit this class and consequently ``self.cell_list`` is a member
of all steppables (please see a chapter on class inheritance from any
Python manual if this looks unfamiliar). Under the hood the
``self.cell_list`` is a handle, or a “pointer”, if you prefer this
terminology, to the C++ object that stores all cells in the simulation.
The content of cell inventory, and cell ordering of cells there is fully
managed by C++ code. We use self.cellList to access C++ cell objects
usually iterating over entire list of cells. The cell in the

.. code-block:: python
for cell in self.cell_list:
.. note::
Functions
****************************

old syntax ``self.cellList`` is still supported
**def __init__(self, frequency=1):** This code runs as the simulation is set up, and, in most cases, you will not need to edit it.

is a pointer to C++ cell object. You can easily see what members C++
cell object has by modifying the step function as follows:
**def start(self):** This is called after cells are created but before the simulation starts, so use it to assign custom cell properties or create `plots <example_plots_histograms.html>`_.

.. code-block:: python
**def step(self, mcs):** Almost everything will happen here. For example, you might grow, divide, or kill your cells here.

def step(self, mcs):
for cell in self.cellList:
print dir(cell)
break
**def on_stop(self):** This runs when you click the stop button.

The result looks as follows:
**def finish(self):** This function is called at the end of the simulation, but it is used very infrequently.
Be careful: it will only run if the simulation reaches the maximum number of steps as specified by the XML attribute ``<Steps>``.

|image5|
**********************************************

*Figure 6 Checking out properties of a cell C++ object*

The dir built-in Python function prints out names of members of any
Python object. Here it printed out members of ``CellG`` class which
represents CC3D cells. We will go over these properties later.

The simplicity of the above code snippets is mainly due to underlying
implementation of ``SteppableBasePy`` class. You can find this class in
``<CC3D\_installation\_dir>/pythonSetupScripts/PySteppables.py``. The
definition of this class goes on for several hundreds lines of code
(clearly a bit too much to present it here). If you are interested in
checking out what members this class has use the dir Python function
again:
A very common line in a Python steppable will read:

.. code-block:: python
def step(self,mcs):
print ('Members of SteppableBasePy class')
print dir(self)
for cell in self.cell_list:
You should know from Python programing manual that ``self`` refers to the
class object. Therefore by printing ``dir(self)`` we are actually printing
Python list of all members of ``cellsortingSteppable`` class. Because
``cellsortingSteppable`` class contains all the functions of ``SteppableBasePy``
class we can inspect this way base class ``SteppableBasePy``. The output of
the above simulation should look as follows:
``cell_list`` is a variable of the ``SteppableBasePy`` class, which
means that every steppable you create will automatically store
every cell it creates in that list.

|image6|
All steppables in CompuCell3D are extensions of SteppableBasePy, so they,
too, can use ``cell_list``.

*Figure 7 Printing all members of SteppableBasePy class*
**If you're new to programming:** The word ``self`` refers to the class we're inside, which would the steppable.
Please see chapters on classes and inheritance from any Python manual if this looks unfamiliar.
Also, if you want a full list of the cell attributes, see `Appendix B <appendix_b.rst>`_.

If you look carefully, you can see that cellList is a member of
``SteppabeBasePy`` class. Alternatively you can study source code of
``SteppablBasePy``.
**If you're a programmer:** Under the hood, the ``self.cell_list`` is a handle, or a “pointer”, to the C++ object that stores all cells in the simulation.
The content of cell inventory and cell ordering of cells there is fully
managed by C++ code.
You can easily see what member variables the C++
cell object has by calling ``dir(cell)`` with one of the cells from the ``self.cell_list`` or by checking out `Appendix B <appendix_b.rst>`_.

One of the goals of this manual is to teach you how to effectively use
features of ``SteppableBasePy`` class to create complex biological
simulations. This class is very powerful and has many constructs which
make coding simple.
.. note::

.. |image5| image:: images/image6.jpeg
:width: 5.98958in
:height: 0.89583in
.. |image6| image:: images/image7.jpeg
:width: 6.00000in
:height: 2.03125in
Old syntax like ``self.cellList`` is still supported.
2 changes: 1 addition & 1 deletion docs/building_a_wall.rst
@@ -1,4 +1,4 @@
Building a wall (it is going to be terrific. Believe me)
Building a wall at the lattice boundary
========================================================

One of the side effects of the Cellular Potts Model occurring when
Expand Down
17 changes: 10 additions & 7 deletions docs/cc3d_python.rst
Expand Up @@ -7,17 +7,20 @@ How to programmatically control every aspect of the simulation

Most new users start learning CompuCell3D by configuring and editing simulations in `Twedit GUI <https://github.com/CompuCell3D/cc3d-twedit5>`_
and then running them via `Player GUI <https://github.com/CompuCell3D/cc3d-player5>`_.
Starting from version 4.4.0 CompuCell3D is shipped with an upgraded Python API that allows modelers to bypass
Starting from version 4.4.0, CompuCell3D is shipped with an upgraded Python API that allows modelers to bypass
both Player and Twedit and configure simulations
- or even an ensemble of multiple simulations running in parallel and informing each other - directly from a
single Python script. This way of constructing simulations gives users full control over every aspect of the simulation.
The material presented in this section explains how to configure entire simulation in pure Python
(we are not using XML), how to write and control main CC3D loop where at every step we call Potts algorithm, how to
run multiple concurrent and interacting simulations, how to specify and control visualization and how to use
Jupyter notebook.
The material presented in this section explains:

In a sense this tutorial shows an alternative way of specifying and running simulations that bypasses legacy tools like
Twedit++ and Player and gives user full programmatic control of the CC3D simulation.
* how to configure the entire simulation in pure Python (*i.e.*, without XML)
* how to write and control the main CC3D loop where, at every step, we call the `Potts algorithm <potts.html>`_
* how to run multiple concurrent and interacting simulations
* how to specify and control visualization
* how to use Jupyter notebook.

In a sense, this tutorial shows an alternative way of specifying and running simulations that bypasses legacy tools like
Twedit++ and Player. It gives you full programmatic control of the CC3D simulation.

We recommend that you use PyCharm or VS Code if you decide to work with the API presented in this section

Expand Down
4 changes: 2 additions & 2 deletions docs/conf.py
Expand Up @@ -55,9 +55,9 @@
# built documents.
#
# The short X.Y version.
version = u'4.4.1'
version = u'4.5.0'
# The full version, including alpha/beta/rc tags.
release = u'4.4.1'
release = u'4.5.0'

# The language for content autogenerated by Sphinx. Refer to documentation
# for a list of supported languages.
Expand Down
133 changes: 120 additions & 13 deletions docs/contact_plugin.rst
@@ -1,29 +1,94 @@
Contact Plugin
--------------

Contact plugin implements computations of adhesion energy between neighboring cells.
Relevant Examples:
* `Cell Sorting <example_cell_sorting.html>`_

**********************************************

Together with volume constraint contact energy is one of the most
commonly used energy terms in the GGH Hamiltonian. In essence it
The contact plugin implements computations of adhesion energy between neighboring cells.

Together with volume constraint, contact energy is one of the most
commonly used plugins. In essence, it
describes how cells "stick" to each other.

The explicit formula for the energy is:
Contact energy is contrived.
It is merely a way to replicate the properties of a cell's membrane, the bindings of the nano-structures on its surface, and its environment (the Medium).

Two cell types that have *high* contact energy will not "want to" adhere to each other.
If possible, those cells may separate, effectively reducing the total energy to stay in that position.
Conversely, *low* contact energy "encourages" cell types to bind.
As contact energy is lowered, it also increases the surface of the contact.

Contact energy is constantly re-calculated each time a cell's surface changes.

We define contact energy as a matrix that compares each cell type against each other cell type.
In CC3DML, try changing the contact energy between CellA and CellB to ``2``.
Keep everything else the same, then click Play |Play| to see the result.

.. code-block:: xml
<Plugin Name="Contact">
<Plugin Name="Contact">
<Energy Type1="CellA" Type2="CellA">10</Energy>
<Energy Type1="CellB" Type2="CellB">10</Energy>
<Energy Type1="CellA" Type2="CellB">2</Energy> <!--Here!-->
<Energy Type1="CellA" Type2="Medium">10</Energy>
<Energy Type1="CellB" Type2="Medium">10</Energy>
<Energy Type1="Medium" Type2="Medium">10</Energy> <!--Medium-Medium energy doesn't really matter-->
<NeighborOrder>2</NeighborOrder>
</Plugin>
This should stimulate cell mixing, the opposite of `cell Sorting <example_cell_sorting.html>`_.


What Does Energy Mean?
-------------------------------------------------

When tweaking your parameters, it's most important to pay attention to their **relative values**.

However, the actual value of each parameter—whether you set ``10`` or ``1``—affects the contact behaviors, too, because of *temperature*.

Each time a cell flexes its surface, a decision is made about which pixels in the lattice will be copied over to new lattice sites.
We discussed this in the `Potts page <potts.html>`_.

The probability of accepting a lattice site copy is:

.. math::
:nowrap:
\begin{eqnarray}
E_{adhesion} = \sum_{i,j,neighbors} J\left ( \tau_{\sigma(i)},\tau_{\sigma(j)} \right )\left ( 1-\delta_{\sigma(i), \sigma(j)} \right )
P = e^{-\bigtriangleup H/T}
\end{eqnarray}
where ``i`` and ``j`` label two neighboring lattice sites :math:`\sigma`'s denote cell
Ids, :math:`\tau`'s denote cell types .
where :math:`\bigtriangleup H = H_f - H_i`.
:math:`H_f` is the new energy after a potential lattice site copy.
:math:`H_i` is the current energy of the system at step *i*.
The new energy, :math:`H_i`, may be higher, lower, or the same after a site copy.
Meanwhile, :math:`T` is the temperature (or fluctuation amplitude *or* membrane fluctuation level *or* noise level).
A high temperature lowers the probability that site copy attempts will be made.


Calculating Contact Energy
-------------------------------------------------

The explicit formula for contact energy is:

.. math::
:nowrap:
\begin{eqnarray}
E_{contact} = \sum_{i,j,neighbors} J\left ( \tau_{\sigma(i)},\tau_{\sigma(j)} \right )\left ( 1-\delta_{\sigma(i), \sigma(j)} \right )
\end{eqnarray}
where ``i`` and ``j`` label two neighboring lattice sites :math:`\sigma`'s denote cell
IDs, and each :math:`\tau` denotes a cell type.
Note that ``i`` and ``j`` are are at separate coordinates; for example, ``i = (0,4)`` and ``j = (1,4)``.

In the above formula, we need to differentiate between cell types and
cell Ids. This formula shows that cell types and cell Ids **are not the
same**. The ``Contact`` plugin in the ``.xml`` file, defines the energy per unit
cell IDs. This formula shows that cell types and cell IDs **are not the
same** and the expression in parentheses involving delta function ensures that we do not included contact surfaces
of voxels that belong to the same cell. The ``Contact`` plugin in the ``.xml`` file, defines the energy per unit
area of contact between cells of different types (:math:`J\left ( \tau_{\sigma(i)},\tau_{\sigma(j)} \right )`) and the interaction
range ``NeighborOrder`` of the contact:

Expand All @@ -37,9 +102,51 @@ range ``NeighborOrder`` of the contact:
</Plugin>
In this case, the interaction range is ``2``, thus only up to second nearest
neighbor pixels of a pixel undergoing a change or closer will be used to calculate
In this case, the interaction range is ``2``. Thus, only up to second-nearest
neighbor pixels of a given pixel undergoing a change will be used to calculate
contact energy change. ``Foam`` cells have contact energy per unit area of ``3``
and ``Foam`` and ``Medium`` as well as Medium and Medium have contact energy of
``0`` per unit area. For more information about contact energy calculations
see “Introduction to CompuCell3D”
``0`` per unit area.

----------------------------------------------------

Suppose we have **8 cells** on the checkboard and **4 sides each**.
If the cells are in a checkerboard pattern, they have the most possible contact energy with the Medium.

.. image:: images/cellsort_legend.png
:alt: A legend for color coding of the two cell types.

.. image:: images/cellsort_checkerboard.png
:alt: An arrangement of 8 cells in a checkboard pattern. The Medium fills in the empty space.

So, the total contact energy is:

.. math::
:nowrap:
\begin{eqnarray}
H_{contact} = 4 \times 8 \times J_{Cell-to-Medium}
\end{eqnarray}
where :math:`J_{Cell-to-Medium}` denotes the contact energy between the base cell type, ``Cell``, and the Medium.

How would you calculate the contact energy if all the cells were touching the medium as little as possible?

.. image:: images/cellsort_compact.png
:alt: An arrangement of 8 cells in a 3x3 area with one corner missing.

.. math::
:nowrap:
\begin{eqnarray}
H_{contact} = 12 \times J_{Cell-to-Medium} + 10 \times J_{Cell-to-Cell}
\end{eqnarray}
Again, count the number of contact surfaces and multiply them by the respective contact energies for those cell types.

.. note::

CompuCell3D never computes total contact energy. Such calculation would be costly and instead CC3D computes a change in Contact energy due to pixel c opy and such calculations are local because we are only examining pixels in a local neighborhood of the change pixel

.. |Play| image:: images/icons/play.png
:height: 14px
39 changes: 39 additions & 0 deletions docs/energy_function_calculator.rst
@@ -0,0 +1,39 @@
How to Output Energy Changes
======================================================

**EnergyFunctionCalculator**: allows you to output statistical data from the simulation for further analysis.

You can define this in XML like so:

.. code-block:: xml
<Potts>
<Dimensions x="256" y="256" z="1"/>
<Steps>100000</Steps>
<Temperature>5</Temperature>
<NeighborOrder>2</NeighborOrder>
<EnergyFunctionCalculator Type="Statistics">
<OutputFileName Frequency="10">statData.txt</OutputFileName>
<OutputCoreFileNameSpinFlips Frequency="1" GatherResults="" OutputAccepted="" OutputRejected="" OutputTotal=""/>
</EnergyFunctionCalculator>
</Potts>
.. note::

CC3D has the option to run in parallel mode, but
output from the energy calculator will only work when running in single-CPU mode.

The ``OutputFileName`` tag is used to specify the name of the file to which
CompuCell3D will write average changes in energies returned by each
plugin along with the corresponding standard deviations.
``Frequency`` controls how often to record this data. In the above example, we would write data every 10 MCS.

Furthermore, the ``OutputCoreFileNameSpinFlips`` tag is used to tell
CompuCell3D to output the energy change for every plugin's pixel-copy operations.

Finally, ``GatherResults=””`` will ensure
that there is only one file written for accepted (``OutputAccepted``),
rejected (``OutputRejected``), and both accepted and rejected (``OutputTotal``) pixel copies.
If you do not specify ``GatherResults``, CompuCell3D will output
separate files for different MCS's, and depending on the Frequency, you
may end up with many files in your directory.

0 comments on commit 6c72b00

Please sign in to comment.