Computing Heterotypic Boundary Length

Heterotypic boundary surface (or length in 2D) is total surface of contact between all cells of two types. For example when you have 2 cells of type 1 and 100 cells of type 2 the heterotypic surface between thw two will be a sum of all contact surfaces between the two types. In this example we are not going to show every step how we generate the steppable using Twedit. We have shown this earlier and here we will concentrate on the actual code.

This example is a bit more advanced but we will explain clearly every line of code.

The module that we generated is called HeterotypicBoundaryLength. We then click Configure and Generate in the CMake Gui and start writing actual code. We will first implement function that walks over entire lattice and computes heterotypic surface (or length in 2D) between all cells of different types.

All C++ files can be found in DeveloperZone/Demos/HeterotypicBoundarySurface and Python bindings are , as usual in DeveloperZone/pyinterface/CompuCellExtraModules/CompuCellExtraModules.i. The simulation example that uses our newly created module is in DeveloperZone/Demos/HeterotypicBoundarySurface

Here is the header file for our new steppable:

#ifndef HETEROTYPICBOUNDARYLENGTHSTEPPABLE_H
#define HETEROTYPICBOUNDARYLENGTHSTEPPABLE_H
#include <CompuCell3D/CC3D.h>
#include "HeterotypicBoundaryLengthDLLSpecifier.h"

namespace CompuCell3D {

  template <class T> class Field3D;
  template <class T> class WatchableField3D;

    class Potts3D;
    class Automaton;
    class BoundaryStrategy;
    class CellInventory;
    class CellG;

  class HETEROTYPICBOUNDARYLENGTH_EXPORT HeterotypicBoundaryLength : public Steppable {

    WatchableField3D<CellG *> *cellFieldG;

    Simulator * sim;
    Potts3D *potts;
    CC3DXMLElement *xmlData;
    Automaton *automaton;
    BoundaryStrategy *boundaryStrategy;
    CellInventory * cellInventoryPtr;
    Dim3D fieldDim;

  public:

    HeterotypicBoundaryLength ();

    virtual ~HeterotypicBoundaryLength ();

    // new methods and members

    std::map<unsigned int, double> typePairHTSurfaceMap;

    unsigned int typePairIndex(unsigned int cellType1, unsigned int cellType2);
    void calculateHeterotypicSurface();
    double getHeterotypicSurface(unsigned int cellType1, unsigned int cellType2);

    // SimObject interface

    virtual void init(Simulator *simulator, CC3DXMLElement *_xmlData=0);

    virtual void extraInit(Simulator *simulator);

    //steppable interface

    virtual void start();

    virtual void step(const unsigned int currentStep);

    virtual void finish() {}

    //SteerableObject interface

    virtual void update(CC3DXMLElement *_xmlData, bool _fullInitFlag=false);

    virtual std::string steerableName();

     virtual std::string toString();


  };

};

#endif

we added few methods and one class member there:

// new methods and members

std::map<unsigned int, double> typePairHTSurfaceMap;

unsigned int typePairIndex(unsigned int cellType1, unsigned int cellType2);
void calculateHeterotypicSurface();
double getHeterotypicSurface(unsigned int cellType1, unsigned int cellType2);

The typePairHTSurfaceMap is a dictionary (map) that will store heterotypic boundary surface between different cell types. Notice that we will encode pair of cell types as a single unsigned integer (hence a key to the dictionary is unsigned integer). To do this we will use convenience function unsigned int typePairIndex(unsigned int cellType1, unsigned int cellType2) that takes as its arguments two unsigned integers that denote cell type 1 and cell type 2. Here is the implementation of this function:

unsigned int HeterotypicBoundaryLength::typePairIndex(unsigned int cellType1, unsigned int cellType2) {
    return 256 * cellType2 + cellType1;
}

we take advantage of the fact that the number of cell types in CC3D is limited to 256 and the index this function returns looks analogous to the index you woudl use to access a matrix if you were to store a matrix as 1D array.

Next we have two functions calculateHeterotypicSurface() that computed actual total heterotypic surface between all cell types and double getHeterotypicSurface(unsigned int cellType1, unsigned int cellType2) that given two types it returns a boundary between them.

Let’s start analyzing code for calculateHeterotypicSurface function:

 1void HeterotypicBoundaryLength::calculateHeterotypicSurface() {
 2
 3    unsigned int maxNeighborIndex = this->boundaryStrategy->getMaxNeighborIndexFromNeighborOrder(1);
 4    Neighbor neighbor;
 5
 6    CellG *nCell = 0;
 7
 8    this->typePairHTSurfaceMap.clear();
 9
10    // note: unit surface is different on a hex lattice. if you are runnign
11    // this steppable on hex lattice you need to adjust it. Remember that on hex lattice unit length and unit surface have different values
12    double unit_surface = 1.0;
13
14    cerr << "Calculating HTBL for all cell type combinations" << endl;
15
16    for (unsigned int x = 0; x < fieldDim.x; ++x)
17        for (unsigned int y = 0; y < fieldDim.y; ++y)
18            for (unsigned int z = 0; z < fieldDim.z; ++z) {
19                Point3D pt(x, y, z);
20                CellG *cell = this->cellFieldG->get(pt);
21
22                unsigned int cell_type = 0;
23                if (cell) {
24                    cell_type = (unsigned int)cell->type;
25                }
26
27                for (unsigned int nIdx = 0; nIdx <= maxNeighborIndex; ++nIdx) {
28                    neighbor = boundaryStrategy->getNeighborDirect(const_cast<Point3D&>(pt), nIdx);
29                    if (!neighbor.distance) {
30                        //if distance is 0 then the neighbor returned is invalid
31                        continue;
32                    }
33
34                    nCell = this->cellFieldG->get(neighbor.pt);
35                    unsigned int n_cell_type = 0;
36                    if (nCell) {
37                        n_cell_type = (unsigned int)nCell->type;
38                    }
39
40                    if (nCell != cell) {
41                        unsigned int pair_index_1 = typePairIndex(cell_type, n_cell_type);
42                        unsigned int pair_index_2 = typePairIndex(n_cell_type, cell_type);
43                        this->typePairHTSurfaceMap[pair_index_1] += unit_surface;
44                        if (pair_index_1 != pair_index_2) {
45                            this->typePairHTSurfaceMap[pair_index_2] += unit_surface;
46                        }
47
48                    }
49
50                }
51
52            }
53
54}

We will be iterating over lattice pixels. Every lattice pixel has neighbors of different order but 1-st order neighbors are simply adjacent pixels. BoundaryStrategy is an object that facilitates iteration over pixel neighbors and it also keeps track of boundary conditions, pixels, adjacent to the boundary etc. so that you can write a simpler code. All we need to do to iterate over 1-st order pixel neighbors is to know what is the maximum number of them and this is what we do in this line:

unsigned int maxNeighborIndex = this->boundaryStrategy->getMaxNeighborIndexFromNeighborOrder(1);

We get maximum index of a 1-st order pixel (BoundaryStrategy keeps them in a vector and we are getting max index of this vector). On 2D. cartesian lattice there could be up to 4 such neighbors hence the max vector index is 3 (we start counting from 0).

We next clear the map where we store our results because each time we call this function it wil be incrementing the values so if we did not clear we would be starting counting from different value that zero.

At line 16 we start triple loop where we iterate over all lattice pixels. This might not be the most efficient method but it is the simplest to code.

In lines 19 and 20 we get a cell that resides at a given pixel. If the cell pointer returned is NULL we are dealing with Medium and cell and this is why in lines 23-25 we check if the cell is different than NULL (if (cell)) before accessing its type. If it is null that we do not execute line 24 and the cell type is 0 as it should be for the Medium.

At line 27 we start iterating over neighbors of the current pixel (Point3D pt(x, y, z)). This is where we do actual calculations. Code in line 28 fetches one of the neighbor of pixel pt(x, y, z). In line 30 we check if this neighbor is a valid one (e.g. if you are at the edge of the lattice we may get pixel that is outside of the lattice and then if neighbor.distance is zero we know we are dealing with invalid pixel), hence in the line 31 we skip the rest of the loop. If, however, the pixel is valid then we get a cell that resides at the neighboring pixel ( line 34):

nCell = this->cellFieldG->get(neighbor.pt);

In lines 35-38 we extract cell type of neighbor cell , again we have to be mindful of Medium as we did in lines 22-25.

Finally, lines 41-45 contain actual code that increments boundary surface between two cell types. This code runs only if nCell and cell i.e. cells belonging to adjacent pixels are different cells. In this case we compute index for type of nCell and type of cell (unsigned int pair_index_1 = typePairIndex(cell_type, n_cell_type);)

and increment appropriate entry in the this->typePairHTSurfaceMap - lines 43. Notice that we also permute cell types in call to typePairIndex - line 44-45. so that when we access boundary length between cell type 1 and 2 it will give us the same value as between cell types 2 and 1. But we do this only when the two types are different

Obviously, we are double-counting and we correct this in the function that returns heterotypic surfaces:

double HeterotypicBoundaryLength::getHeterotypicSurface(unsigned int cellType1, unsigned int cellType2) {
    unsigned int pair_index = typePairIndex(cellType1, cellType2);

    double heterotypic_surface = this->typePairHTSurfaceMap[pair_index]/2.0;

    return heterotypic_surface;
}

Running the Simulation with Heterotypic Surface Calculator

The simulation code is quite easy to write as it follows the same pattern that we encountered in the previous chapter where we introduced Python bindings to the C++ steppable. We start with an XML file:

 1<CompuCell3D Revision="20190604" Version="4.0.0">
 2   <Potts>
 3      <!-- Basic properties of CPM (GGH) algorithm -->
 4      <Dimensions x="256" y="256" z="1"/>
 5      <Steps>100000</Steps>
 6      <Temperature>10.0</Temperature>
 7      <NeighborOrder>1</NeighborOrder>
 8   </Potts>
 9
10   <Plugin Name="CellType">
11      <!-- Listing all cell types in the simulation -->
12      <CellType TypeId="0" TypeName="Medium"/>
13      <CellType TypeId="1" TypeName="A"/>
14      <CellType TypeId="2" TypeName="B"/>
15   </Plugin>
16
17   <Plugin Name="Volume">
18      <VolumeEnergyParameters CellType="A" LambdaVolume="2.0" TargetVolume="50"/>
19      <VolumeEnergyParameters CellType="B" LambdaVolume="2.0" TargetVolume="50"/>
20   </Plugin>
21
22   <Plugin Name="CenterOfMass">
23      <!-- Module tracking center of mass of each cell -->
24   </Plugin>
25
26   <Plugin Name="Contact">
27      <!-- Specification of adhesion energies -->
28      <Energy Type1="Medium" Type2="Medium">10.0</Energy>
29      <Energy Type1="Medium" Type2="A">10.0</Energy>
30      <Energy Type1="Medium" Type2="B">10.0</Energy>
31      <Energy Type1="A" Type2="A">10.0</Energy>
32      <Energy Type1="A" Type2="B">10.0</Energy>
33      <Energy Type1="B" Type2="B">10.0</Energy>
34      <NeighborOrder>4</NeighborOrder>
35   </Plugin>
36
37   <Steppable Type="UniformInitializer">
38      <!-- Initial layout of cells in the form of rectangular slab -->
39      <Region>
40         <BoxMin x="51" y="51" z="0"/>
41         <BoxMax x="204" y="204" z="1"/>
42         <Gap>0</Gap>
43         <Width>7</Width>
44         <Types>A,B</Types>
45      </Region>
46   </Steppable>
47
48   <Steppable Type="HeterotypicBoundaryLength"/>
49
50</CompuCell3D>

It is Twedit-generated XML file that has basic energy terms (Volume and Contact Constraints) plus initializer and at the end in line 48 we add our new HeterotypicBoundaryLength steppable. Notice that this is a one-line call because we are not really passing any parameters to the steppable from the XML and our update(CC3DXMLElement *_xmlData, bool _fullInitFlag=false) method does not contain any code that parses XML.

Note

It is important that every module (steppable, plugin) that you develop in C++ be instantiated in XML. Otherwise it will not be loaded and you will not be able to use it from Python. You can, however, write Python code that will properly load and initialize your module but this approach is way more complex than adding a simple line or lines in the XML.

In our example even if we add <Steppable Type="HeterotypicBoundaryLength"/> to the XML we will not see any calculations being done. Why? Because start and step functions are empty:

void HeterotypicBoundaryLength::start(){
}

void HeterotypicBoundaryLength::step(const unsigned int currentStep){
}

void HeterotypicBoundaryLength::update(CC3DXMLElement *_xmlData, bool _fullInitFlag){

    //PARSE XML IN THIS FUNCTION

    //For more information on XML parser function please see CC3D code or lookup XML utils API

    automaton = potts->getAutomaton();

    ASSERT_OR_THROW("CELL TYPE PLUGIN WAS NOT PROPERLY INITIALIZED YET. MAKE SURE THIS IS THE FIRST PLUGIN THAT YOU SET", automaton)

    //boundaryStrategy has information about pixel neighbors
    boundaryStrategy=BoundaryStrategy::getInstance();

}

We left those implementations empty on purpose. We wanted to show you how you can use steppable to implement functionality that gets called on-demand from Python code. Let us now look at the Python code:

 1from cc3d.core.PySteppables import *
 2from cc3d.cpp import CompuCellExtraModules
 3
 4
 5class HeterotypicBoundarySurfaceSteppable(SteppableBasePy):
 6
 7    def __init__(self, frequency=1):
 8        SteppableBasePy.__init__(self, frequency)
 9        self.htbl_steppable_cpp = None
10
11    def start(self):
12        self.htbl_steppable_cpp = CompuCellExtraModules.getHeterotypicBoundaryLength()
13
14    def step(self, mcs):
15        self.htbl_steppable_cpp.calculateHeterotypicSurface()
16
17        print(' HTBL between type 1 and 2 is ',
18              self.htbl_steppable_cpp.getHeterotypicSurface(1, 2))
19
20        print(' HTBL between type 2 and 1 is ',
21              self.htbl_steppable_cpp.getHeterotypicSurface(1, 2))
22
23        print(' HTBL between type 1 and 1 is ',
24              self.htbl_steppable_cpp.getHeterotypicSurface(1, 1))
25
26        print(' HTBL between type 0 and 1 is ',
27              self.htbl_steppable_cpp.getHeterotypicSurface(0, 1))
28
29        print('THIS ENTRY DOES NOT EXIST. HTBL between type 3 and 20 is ',
30              self.htbl_steppable_cpp.getHeterotypicSurface(3, 20))

At the top we import CompuCellExtraModules that SWIG generates. This Python contains contains Python-wrapper for our HeterotypicBoundaryLength C++ steppable. In line 12 we fetch a reference to the steppable and store it in class-accessible self.htbl_steppable_cpp variable. In Python steppable step function we call C++ function that calculates heterotypic surface (line 15). The next series of print statements fetches results fo the calculations - see lines 18 21, etc… . The output looks as follows:

htbl_output

Note that heterotypic surface between types 1 and 2 is the same as between 2 and 1. This is why, earlier in the C++ code we included two pair indexes:

unsigned int pair_index_1 = typePairIndex(cell_type, n_cell_type);
unsigned int pair_index_2 = typePairIndex(n_cell_type, cell_type);

Notice also that if we try to access heterotypic surface between types 3 and 20 (none of those types exist in our simulation) we get back 0.0. Why? The answer has to do with the behavior of C++ std::map container. If we try to use a key that does not exist in the map the C++ will add this key and initialize the value of the key value-pair to the whatever default constructor of the value type is. In our case our map container has the following type : std::map<unsigned int, double> typePairHTSurfaceMap so that key is of unsigned int type and value is of double type. Any modern C++ compiler will put 0.0 as a default value for objects of type double. If you are familiar with Python defaultdict class that is a member of standard collections package than you can see similarities. C++ std::map behaves in similar way to the defaultdict Therefore when we access typePairHTSurfaceMap using key that does not exist C++ will insert this key into typePairHTSurfaceMap and set value to 0.0. This is also a reason why the following code works at all:

this->typePairHTSurfaceMap[pair_index_1] += unit_surface

This brings us to the last remark we want to make regarding the C++ code. Why are we using unit_surface and not 1.0? It has to do with various lattices that CC3D supports. For Cartesian 2D or 3D lattice unit surface has value 1.0. However, hexagonal lattices are constructed in such a way that the volume of a voxel is constrained to be 1.0. Therefore from geometry constraints it follows that in 2D on hex lattice unit surface (or length) is sqrt(2.0 / (3.0*sqrt(3.0))) which is approx equal to 0.6204 and in 3D it is 8.0 / 12.0*sqrt(2.0)*pow(9.0 / (16.0*sqrt(3.0)), 1.0 / 3.0)*pow(9.0 / (16.0*sqrt(3.0)), 1.0 / 3.0) which is approx equal to 0.445

For more information please see http://www.compucell3d.org/BinDoc/cc3d_binaries/Manuals/HexagonalLattice.pdf as well as in the code of the BoundaryStrategy class method LatticeMultiplicativeFactors BoundaryStrategy::generateLatticeMultiplicativeFactors(LatticeType _latticeType, Dim3D dim) in CompuCell3D/core/CompuCell3D/Boundary/BoundaryStrategy.cpp

For completeness we also show SWIG file that was used to generate the wrapper :

%module ("threads"=1) CompuCellExtraModules

%include "typemaps.i"

%include <windows.i>

%{

#include "ParseData.h"
#include "ParserStorage.h"
#include <CompuCell3D/Simulator.h>
#include <CompuCell3D/Potts3D/Potts3D.h>

#include <BasicUtils/BasicClassAccessor.h>
#include <BasicUtils/BasicClassGroup.h> //had to include it to avoid problems with template instantiation


// ********************************************* PUT YOUR PLUGIN PARSE DATA AND PLUGIN FILES HERE *************************************************

#include <SimpleVolume/SimpleVolumePlugin.h>
#include <VolumeMean/VolumeMean.h>

//AutogeneratedModules1 - DO NOT REMOVE THIS LINE IT IS USED BY TWEDIT TO LOCATE CODE INSERTION POINT
//HeterotypicBoundaryLength_autogenerated

#include <HeterotypicBoundaryLength/HeterotypicBoundaryLength.h>

//GrowthSteppable_autogenerated

#include <GrowthSteppable/GrowthSteppable.h>


// ********************************************* END OF SECTION ********************************** ************************************************


//have to include all  export definitions for modules which are arapped to avoid problems with interpreting by swig win32 specific c++ extensions...
#define SIMPLEVOLUME_EXPORT
#define VOLUMEMEAN_EXPORT
//AutogeneratedModules2 - DO NOT REMOVE THIS LINE IT IS USED BY TWEDIT TO LOCATE CODE INSERTION POINT
//HeterotypicBoundaryLength_autogenerated
#define HETEROTYPICBOUNDARYLENGTH_EXPORT
//GrowthSteppable_autogenerated
#define GROWTHSTEPPABLE_EXPORT

#include <iostream>

using namespace std;
using namespace CompuCell3D;

%}

// C++ std::string handling
%include "std_string.i"

// C++ std::map handling
%include "std_map.i"

// C++ std::map handling
%include "std_set.i"

// C++ std::vector handling
%include "std_vector.i"

//have to include all  export definitions for modules which are arapped to avoid problems with interpreting by swig win32 specific c++ extensions...
#define SIMPLEVOLUME_EXPORT
#define VOLUMEMEAN_EXPORT

//AutogeneratedModules3 - DO NOT REMOVE THIS LINE IT IS USED BY TWEDIT TO LOCATE CODE INSERTION POINT
//HeterotypicBoundaryLength_autogenerated
#define HETEROTYPICBOUNDARYLENGTH_EXPORT
//GrowthSteppable_autogenerated
#define GROWTHSTEPPABLE_EXPORT

%include <BasicUtils/BasicClassAccessor.h>
%include <BasicUtils/BasicClassGroup.h> //had to include it to avoid problems with template instantiation
%include "ParseData.h"
%include "ParserStorage.h"


// ********************************************* PUT YOUR PLUGIN PARSE DATA AND PLUGIN FILES HERE *************************************************
// REMEMBER TO CHANGE #include to %include
%include <SimpleVolume/SimpleVolumePlugin.h>
// %include <SimpleVolume/SimpleVolumeParseData.h>
// THIS IS VERY IMORTANT STETEMENT WITHOUT IT SWIG will produce incorrect wrapper code which will compile but will not work
using namespace CompuCell3D;

%inline %{
   SimpleVolumePlugin * reinterpretSimpleVolumePlugin(Plugin * _plugin){
      return (SimpleVolumePlugin *)_plugin;
   }

   SimpleVolumePlugin * getSimpleVolumePlugin(){
         return (SimpleVolumePlugin *)Simulator::pluginManager.get("SimpleVolume");
    }

%}

%include <VolumeMean/VolumeMean.h>

%inline %{
   VolumeMean * reinterpretVolumeMean(Steppable * _steppable){
      return (VolumeMean *)_steppable;
   }

   VolumeMean * getVolumeMeanSteppable(){
         return (VolumeMean *)Simulator::steppableManager.get("VolumeMean");
    }
%}

//AutogeneratedModules4 - DO NOT REMOVE THIS LINE IT IS USED BY TWEDIT TO LOCATE CODE INSERTION POINT
//HeterotypicBoundaryLength_autogenerated
%include <HeterotypicBoundaryLength/HeterotypicBoundaryLength.h>
%inline %{

 HeterotypicBoundaryLength * getHeterotypicBoundaryLength(){
      return (HeterotypicBoundaryLength *)Simulator::steppableManager.get("HeterotypicBoundaryLength");
   }
%}

//GrowthSteppable_autogenerated
%include <GrowthSteppable/GrowthSteppable.h>
%inline %{

 GrowthSteppable * getGrowthSteppable(){
      return (GrowthSteppable *)Simulator::steppableManager.get("GrowthSteppable");
   }
%}

// ********************************************* END OF SECTION ********************************** ************************************************