# Appendix¶

This section covers some of the computational algorithms we use in CC3D

- Calculating Inertia tensor In CC3D
- Calculating shape constraint of a cell – elongation term
- Forward Euler method for solving PDE’s in CompuCell3D.
- Calculating center of mass when using periodic boundary conditions.
- Dividing cluster cells
- Command line options of CompuCell3D
- Managing CompuCell3D simulations (CC3D project files)

## Calculating Inertia tensor In CC3D¶

For each cell the inertia tensor is defined as follows:

where index \(i\) denotes `i-th`

pixel of a given cell and \(x_i\),
\(y_i\) and \(z_i\) are coordinates of that pixel in a given
coordinate frame.

**Figure 3:** Cell and coordinate system passing through center of mass of a cell.
Notice that as cell changes shape the position of center of mass moves.

**Figure 4:** Cell and its coordinate frame in which we calculate inertia tensor

In Figure 4 we show one possible coordinate frame in which one can calculate inertia tensor. If the coordinate frame is fixed calculating components of inertia tensor for cell gaining or losing one pixel is quite easy. We will be adding and subtracting terms like \(y_i^2+z_i^2\) or \(x_i z_i\).

However, in CompuCell3D we are mostly interested in knowing tensor of
inertia of a cell with respect to `xyz`

coordinate frame with origin at
the center of mass (*COM*) of a given cell as shown in Figure 3. Now, to
calculate such tensor we cannot simply add or subtract terms \(y_i^2+z_i^2\) or \(x_i z_i\) to
account for lost or gained pixel. If a cell gains or loses a pixel its
`COM`

coordinates change. If so then all the \(x_i\),:math:y_i, \(z_i\)
coordinates that appear in the inertia tensor
expression will have different values. Thus, for each change in cell shape
(gain or loss of pixel) we would have to recalculate inertia tensor from
scratch. This would be quite time consuming and would require us to keep
track of all the pixels belonging to a given cell. It turns out, however,
that there is a better way of keeping track of inertia tensor for cells.
We will be using parallel axis theorem to do the calculations. Parallel
axis theorem states that if \(I_{COM}\) is a moment of inertia with
respect to axis passing through center of mass then we can calculate
moment of inertia with respect to any parallel axis to the one passing
through the COM by using the following formula:

where \(I_{xx}\) denotes moment of inertia with respect to `x`

axis passing through
center of mass, \(I_{x'x'}\) is a moment of inertia with respect to some other axis parallel to
the `x`

, `d`

is the distance between the axes two and `M`

is mass of the cell.

Let us now draw a picture of a cell gaining one pixel:

**Figure 5:** Cell gaining one pixel.d denotes a distance from origin of a fixed frame
of reference to a center of mass of a cell before cell gains new pixel.
\(d_{new}\) denotes same distance but after cell gains new pixel

Now using parallel axis theorem we can write expression for the moment of inertia after cell gains one pixel the following that:

where, as before, \(I_{xx}^{new}\) denotes moment of inertia of a cell with new pixel with
respect to `x`

axis passing through center of mass, \(I_{x'x'}^{new}\) is a moment of
inertia with respect to axis parallel to the `x`

axis passing through
center of mass, \(d_{new}\) is the distance between the axes and
\(V+1\) is volume of the cell **after** it gained one pixel. Now let us
rewrite above equation by adding ad subtracting \(Vd^2\) term:

Therefore we have found an expression for moment of inertia passing
through the center of mass of the cell with additional pixel. Note that
this expression involves moment of inertia but for the old cell (*i.e*.
the original cell, not the one with extra pixel). When we add new pixel
we know its coordinates and we can also easily calculate \(d_new\) .
Thus, when we need to calculate the moment of inertia for new cell
instead of performing summation as given in the definition of the
inertia tensor we can use much simpler expression.

This was diagonal term of the inertia tensor. What about off-diagonal terms? Let us write explicitly expression for \(I_{xy}\) :

where \(x_{COM}\), \(y_{COM}\) denote `x`

and `y`

center of mass positions of the cell,
\(V\) denotes cell volume. In the above formula we have used the fact that:

and similarly for the `y`

coordinate.

Now, for the new cell with additional pixel we have the following relation:

where we have added and subtracted \(V x_{COM}y_{COM}\) to be able to form \(I_{xy}^{old}- \sum_i^{N} x_i y_i+ V x_{COM}y_{COM}\)
on the right hand side of the expression for \(I_{xy}^{new}\) . As it was the case for diagonal element,
calculating off-diagonal of the inertia tensor involves ::math`I_{xy}^{old}` and positions of
center of mass of the cell before and after gaining new pixel. All those
quantities are either known *a priori* (::math`I_{xy}^{old}`) or can be easily calculated
(center of mass position after gaining one pixel).

Therefore, we have shown how we can calculate tensor of inertia for a given cell with respect to a coordinate frame with origin at cell’s center of mass, without evaluating full sums. Such “local” calculations greatly speed up simulations

## Calculating shape constraint of a cell – elongation term¶

The shape of single cell immersed in medium and not subject to too drastic surface or surface constraints will be spherical (circular in 2D). However in certain situation we may want to use cells which are elongated along one of their body axes. To facilitate this we can place constraint on principal lengths of cell. In 2D it is sufficient to constrain one of the principal lengths of cell how ever in 3D we need to constrain 2 out of 3 principal lengths. Our first task is to diagonalize inertia tensor (i.e. find a coordinate frame transformation which brings inertia tensor to a diagonal form)

### Diagonalizing inertia tensor¶

We will consider here more difficult 3D case. The 2D case is described
in detail in M.Zajac, G.L.jones, J,A,Glazier “*Simulating convergent
extension by way of anisotropic differential adhesion”* Journal of
Theoretical Biology **222** (2003) 247–259.

In order to diagonalize inertia tensor we need to solve eigenvalue equation:

or in full form

The eigenvalue equation will be in the form of 3^{rd} order
polynomial. The roots of it are guaranteed to be real. The polynomial
itself can be found either by explicit derivation, using symbolic
calculation or simply in Wikipedia ( http://en.wikipedia.org/wiki/Eigenvalue_algorithm )

so in our case the eigenvalue equation takes the form:

This equation can be solved analytically, again we may use Wikipedia ( http://en.wikipedia.org/wiki/Cubic_function )

Now, the eigenvalues found that way are principal moments of inertia of a cell. That is they are components of inertia tensor in a coordinate frame rotated in such a way that off-diagonal elements of inertia tensor are 0:

In our cell shape constraint we will want to obtain ellipsoidal cells. Therefore the target tensor of inertia for the cell should be tensor if inertia for ellipsoid:

where `a`

, `b`

, `c`

are parameters describing the surface of an ellipsoid:

In other words `a`

, `b`

, `c`

are half lengths of principal axes (they are
analogues of circle’s radius).

Now we can determine semi axes lengths in terms of principal moments of inertia by inverting the following set of equations:

Once we have calculated semiaxes lengths in terms of moments of inertia we can plug –in actual numbers for moment of inertia (the ones for actual cell) and obtain lengths of semiexes. Next we apply quadratic constraint on largest (semimajor) and smallest (seminimor axes). This is what elongation plugin does.

## Forward Euler method for solving PDE’s in CompuCell3D.¶

Note

We present more complete derivations of explicit finite difference scheme for diffusion solver in “Introduction to Hexagonal Lattices in CompuCell3D” (http://www.compucell3d.org/BinDoc/cc3d_binaries/Manuals/HexagonalLattice.pdf).

In CompuCell3D most of the solvers uses explicit schemes (Forward Euler method) to obtain PDE solutions. Thus for the diffusion equation we have:

In a discretetized form we may write:

where to save space we used shorthand notation:

and similarly for other coordinates.

After rearranging terms we get the following expression:

where the sum over index \(i\) goes over neighbors of point \((x,y,z)\) and the neighbors will have the following concentrations: \(c(x+\delta x, t)\), \(c(y+\delta y, t)\), \(c(z+\delta z, t)\) .

## Calculating center of mass when using periodic boundary conditions.¶

When you are running calculation with periodic boundary condition you may end up with situation like in the figure below:

**Figure 6** A connected cell in the lattice edge area – periodic boundary conditions are applied

Clearly, what happens is that simply connected cell is wrapped around the lattice edge so part of it is in the region of high values of x coordinate and the other is in the region where x coordinates have low values. Consequently, a naïve calculation of center of mass position according to:

or in vector form:

would result in being somewhere in the middle of the lattice and obviously outside the cell. A better procedure could be as follows:

Before calculating center of mass when new pixel is added or lost we “shift” a cell and new pixel (gained or lost )to the middle of the lattice do calculations “in the middle of the lattice” and shift back. Now if after shifting back it turns out that center of mass of a cell lies outside lattice position it in the center of mass by applygin a shift equal to the length of the lattice and whose direction should be such that the center of mass of the cell ends up inside the lattice (there is only one such shift and it might be be equal to zero vector).

This is how we do it using mathematical formulas:

First we define shift vector \(\vec{s}\) as a vector difference between vector pointing to center of mass of the cell \(\vec{r}_{COM}\) and vector pointing to (approximately) the middle of the lattice \(\vec{c}\).

Next we shift cell to the middle of the lattice using :

where \(\vec{r'}_{COM}\) denotes center of mass position of a cell after shifting but before adding or subtracting a pixel.

Next we take into account the new pixel (either gained or lost) and calculate center of mass position (for the shifted cell):

Above we have assumed that we are adding one pixel.

Now all that we need to do is to shift back \(\vec{r'}_{COM}^{new}\) by same vector \(\vec{s}\) that brought cell to (approximately) center of the lattice:

We are almost done. We still have to check if \(\vec{r'}_{COM}^{new}\) is inside the lattice. If this is not the case we need to shift it back to the lattice but now we are allowed to use only a vector \(\vec{P}\) whose components are multiples of lattice dimensions (and we can safely restrict to +1 and -1 multiples of the lattice dimensions) . For example we may have:

where \(\vec{x}_{max}\), \(\vec{y}_{max}\), \(\vec{z}_{max}\) are dimensions of the lattice.

There is no cheating here. In the lattice with periodic boundary conditions you are allowed to shift point coordinates a vector whose components are multiples of lattice dimensions.

All we need to do is to examine new center of mass position and form suitable vector \(\vec{P}\).

## Dividing cluster cells¶

While dividing non-clustered cells is straightforward, doing the same for clustered cells is more challenging. To divide non-cluster cell using directional mitosis algorithm we construct a line or a plane passing through center of mass of a cell and pixels of the cell (we are using PixelTracker plugin with mitosis) on one side of the line/plane end up in child cell and the rest stays in parent cell. The orientation of the line/plane can be either specified by the user or we can use CC3D built-in feature to calculate calculate orientation tion of principal axes and divide either along minor or major axis.

With compartmental cells, things get more complicated because: 1) Compartmental cells are composed of many subcells. 2) There can be different topologies of clusters. Some clusters may look “snake-like” and some might be compactly packed blobs of subcells. The algorithm which we implemented in CC3D works in the following way:

- We first construct a set of pixels containing every pixel belonging to a cluster cell. You may think of it as of a single “regular” cell.
- We store volumes of compartments so that we know how big compartments should be after mitosis (they will be half of original volume)
- We calculate center of mass of entire cluster and calculate vector offsets between center of mass of a cluster and center of mass of particular compartments as on the figure below:

**Figure 7**.Vectors \(\vec{o}_1\) and \(\vec{o}_2\) show offsets between center of mass of a
cluster and center of mass particular compartments.

4) We pick division line/plane and for parents and child cells with offsets between cluster center of mass (after mitosis) and center of masses of clusters. We do it according to the formula:

where \(\vec{p}\) denotes offset after mitosis from center of mass of child (parent) clusters, \(\vec{o}\) is orientation vector before mitosis (see picture above) and \(\vec{n}\) is a normalized vector perpendicular to division line/plane. If we try to divide the cluster along dashed line as on the picture below:

**Figure 8**. Division of cell along dashed line. Notice the orientation of \(\vec{n}\) .
The offsets after the mitosis for child and parent cell will be
\(\vec{p}_1=\frac{1}{2}\vec{o}_1\) and \(\vec{p}_2=\vec{o}_2\) as
expected because both parent and child cells will retain their heights
but after mitosis will become twice narrower (cell with grey outer

compartments is a parent cell):

**Figure 9**.Child and parent (the one with grey outer compartments)
cells after mitosis.

The formula given above is heuristic. It gives fairly simple way of assigning pixels of child/parent clusters to cellular compartments. It is not perfect but the idea is to get approximate shape of the cell after the mitosis and as simulation runs cell shape will readjust based on constraints such as adhesion of focal point plasticity. Before continuing with the mitosis we check if center of masses of compartments belong to child/parent clusters. If the center of masses are outside their target pixels we abandon mitosis and wait for readjustment of cell shape at which point mitosis algorithm will pass this sanity check. For certain “exotic” shapes of cluster shapes presented mitosis algorithm may not work well or at all. In this case we would have to write specialized mitosis algorithm.

5) We divide clusters and knowing offsets from child/parent cluster center of mass we assign pixels to particular compartments. The assignment is based on the distance of particular pixel to center of masses of clusters. Pixel is assigned to particular compartment if its distance to the center of mass of the compartment is the smallest as compared to distances between centroids of other compartments. If given compartment has reached its target volume and other compartments are underpopulated we would assign pixels to other compartments based on the closest distance criterion. Although this method may result in some deviation from perfect 50-50 division of compartment volume in most cases after few MCS cells will readjust their volume.

**Figure 10.** CC3D example of compartmental cell division. See also
Demos/CompuCellPythonTutorial/clusterMitosis.

## Command line options of CompuCell3D¶

Although most users run CC3D using Player GUI sometimes it is very convenient to run CC3D using command line options. CC3D allows to invoke Player directly from command line which is convenient because if saves several clicks and if you run many simulations this might be quite convenient.

**Remark:** On Windows we use .bat extension for run scripts and on
Linux/OSX it is .sh. Otherwise all the material in this section applies
to all the platforms.

### CompuCell3D Player Command Line Options¶

The command line options for running simulation with the player are as follows:

```
compucell3d.bat [options]
```

Options are:

`-i <simulation file>`

- users specify `.cc3d`

simulation file they want to run.

`-s <screenshotDescriptionFileName>`

- name of the file containing
description of screenshots to be taken with the simulation. Usually this
file is prepared using Player by switching to different views, clickin
camera button and saving screenshot description file from the Player
File menu.

`-o <customScreenshotDirectoryName>`

- allows users to specify where
screenshots will be written. Overrides default settings.

`--noOutput`

- instructs CC3D not to store any screenshots. Overrides
Player settings.

`--exitWhenDone`

- instructs CC3D to exit at the end of simulation.
Overrides Player settings.

`-h`

, `--help`

- prints command line usage on the screen

Example command may look like (windows):

```
compucell3d.bat –i Demos\Models\cellsort\cellsort_2D\cellsort_2d.cc3d --noOutput
```

or on linux:

```
compucell3d.sh –i Demos/Models/cellsort/cellsort_2D/cellsort_2d.cc3d --noOutput
```

and OSX:

```
compucell3d.command –i Demos/Models/cellsort/cellsort_2D/cellsort_2d.cc3d --noOutput
```

### Running CompuCell3D in a GUI-Less Mode - Command Line Options.¶

Sometimes when you want to run CC3D on a cluster you will have to use runScript.bat which allows running CC3D simulations without invoking GUI. However, all the screenshots will be still stored.

**Remark:** current version of this script does not handle properly
relative paths so it has to be run from the installation directory of
CC3D i.e. you have to cd into this directory prior to runnit
runScript.bat. Another solution is to use full paths.

The output of this script is in the form of vtk files which can be subsequently replayed in the Player (and one can take screenshots then). By default all fields present in the simulation are stored in the vtk file. If users want to remove some of the fields from being stored in the vtk format they have to pass this information in the Python script:

```
CompuCellSetup.doNotOutputField(_fieldName)
```

Storing entire fields (as opposed to storing screenshots) preserves
exact snapshots of the simulation and allows result postprocessing. In
addition to the vtk files runScript stores lattice description file with
`.dml``

extension which users open in the Player (`File->Open Lattice Description Summary File…`

)
if they want to reply generated vtk files.

The format of the command (windows):

```
runScript.bat [options]
```

linux:

```
runScript.sh [options]
```

OSX:

```
runScript.command [options]
```

The command line options for runScript.bat are as follows:

`-i <simulation file>`

- users specify `.cc3d`

simulation file they want to run.

`-c <outputFileCoreName>`

- allows users to specify core name for the vtk
files. The default name for vtk files is `Step`

`-o <customVtkDirectoryName>`

- allows users to specify where vtk files
and the .dml file will be written. Overrides default settings

`-f <frequency>`

or `–outputFrequency=<frequency>`

- allows to specify how
often vtk files are stored to the disk. Those files tend to be quite
large for bigger simulations so storing them every single MCS (default
setting) slows down simulation considerably and also uses a lot of disk
space.

`--noOutput`

- instructs CC3D not to store any output. This option makes
little sense in most cases.

`-h`

, `--help`

- prints command line usage on the screen

Example command may look as follows(windows):

```
runScript.bat –i Demos\CompuCellPythonTutorial\InfoPrinter\cellsort_2D_info_printer.cc3d –f 10 –o Demos\CompuCellPythonTutorial\InfoPrinter\screenshots –c infoPrinter
```

linux:

```
runScript.sh –i Demos/CompuCellPythonTutorial/InfoPrinter/cellsort_2D_info_printer.cc3d –f 10 –o Demos/CompuCellPythonTutorial/InfoPrinter/screenshots –c infoPrinter
```

osx:

```
runScript.command –i Demos/CompuCellPythonTutorial/InfoPrinter/cellsort_2D_info_printer.cc3d –f 10 –o Demos/CompuCellPythonTutorial/InfoPrinter/screenshots –c infoPrinter
```

## Managing CompuCell3D simulations (CC3D project files)¶

Until version `3.6.0`

CompuCell3D simulations were stored as a combination
of Python, CC3DML (XML), and PIF files. Starting with version
`3.6.0`

we introduced new way of managing CC3D simulations by enforcing
that a single CC3D simulation is stored in a folder containing `.cc3d`

project file describing simulation resources (.cc3d is in fact XML),
such as CC3DML configuration file, Python scripts, PIF files,
Concentration files *etc…** and a directory called `Simulation`

where all
the resources reside. The structure of the new-style CC3D simulation is
presented in the diagram below:

**->CellsortDemo**

CellsortDemo.cc3d

->SimulationCellsort.xml

Cellsort.py

CellsortSteppables.py

Cellsort.piff

FGF.txt

Bold fonts denote folders. The benefit of using CC3D project files instead of loosely related files are as follows:

- Previously users had to guess which file needs to be open in CC3D –
CC3DML or Python. While in a well written simulation one can link the
files together in a way that when user opens either one the
simulation would work but, nevertheless, such approach was clumsy and
unreliable. Starting with
`3.6.0`

users open`.cc3d`

file and they don’t have to stress out that CompUCell3D will complain with error message.

Warning

The only way to load simulation in CompuCell3D is to use `.cc3d`

project. We no longer
support previous ways of opening simulations

- All the files specified in the .cc3d project files are copied to the result output directory along with simulation results (uncles you explicitly specify otherwise). Thus, when you run multiple simulations each one with different parameters, the copies of all CC3DML and Python files are stored eliminating guessing which parameters were associated with particular simulations.
- All file paths appearing in the simulation files are relative paths
with respect to main simulation folder. This makes simulations
portable because all simulation resources are contained withing
single folder. In the example above when referring to
`Cellsort.piff`

file from`Cellsort.xml`

you use`Simulation/Cellsort.piff`

. This makes simulations easily exchangeable between collaborators - New style of storing CC3D simulations has also another advantage – it makes graphical management of simulation content and simulation generation very easy. As a matetr of fact new component of CC3D suite – Twedit++ - CC3D edition has a graphical tool that allows for easy project file management and it also has new simulation wizard which allows users to build template of CC3D simulation within less than a minute.

Let’s now look in detail at the structure of `.cc3d`

files (we are using XML syntax here):

```
<Simulation version="3.6.0">
<XMLScript>Simulation/Cellsort.xml</XMLScript>
<PythonScript>Simulation/Cellsort.py</PythonScript>
<Resource Type="Python">Simulation/CellsortSteppables.py</Resource>
<PIFFile>Simulation/Cellsort.piff</PIFFile>
<Resource Type="Field" Copy="No">Simulation/FGF.txt</Resource>
</Simulation>
```

As you can see the structure of the file is quite flat. All that we are
storing there is names of files that are used in the simulation. Two
files have special tags `<XMLFile>`

which specifies name of the CC3DML
file storing “CC3DML portion” of the simulation and `<PythonScript>`

which
specifies main Python script. We have also `<PIFFile>`

tag which is used to
designate PIF files. All other files used in the simulation are referred
to as `Resources`

. For example Python steppable file is a resource of type
“Python” - `<Resource Type="Python">Simulation/CellsortSteppables.py</Resource>`

.
FGF.txt is a resource of type “Field” - `<Resource Type="Field" Copy="No">Simulation/FGF.txt</Resource>`

.
Notice that all the files are specified using paths relative to main simulation directory i.e. w.r.t to the dir in which
`.cc3d`

file resides

As we mentioned before, when you run `.cc3d`

simulation all the files
listed in the project file are copied to result folder. If for
some reason you want to avoid coping of some of the files, simply add
`Copy="No"`

attribute in the tag with file name specification for example”

```
<PIFFile Copy="No">Simulation/Cellsort.piff</PIFFile>
<Resource Type="Field" Copy="No">Simulation/FGF.txt</Resource>
```