On this page This tutorial was generated from the file projects/PottsCrypt2015/test/TestPottsCryptSweepsLiteratePaper.hpp at revision r24335.
Note that the code is given in full at the bottom of the page.

Multiple CPM simulations of a healty crypt with varying dynamic parameters# Introduction# In this test we run multiple simulations if a Cellular Potts model of a
colon crypt. We do this to see the dependance of the simulation
in the dynamic parameters of the CPM

Full details of the computational model can be found in
Osborne (2015) “A Multiscale Model of Colorectal Cancer Using the Cellular Potts Framework”.

This class was used to produce the data for Figure 4.

We begin by including the necessary header files which are the same as for a single simulation.

#include <cxxtest/TestSuite.h>
#include "CellBasedSimulationArchiver.hpp"
#include "TransitCellProliferativeType.hpp"
#include "PottsMeshGenerator.hpp"
#include "CellsGenerator.hpp"
#include "WntConcentration.hpp"
#include "SimpleWntCellCycleModel.hpp"
#include "PottsBasedCellPopulation.hpp"
#include "VolumeConstraintPottsUpdateRule.hpp"
#include "AdhesionPottsUpdateRule.hpp"
#include "SloughingCellKiller.hpp"
#include "OnLatticeSimulation.hpp"
#include "CellProliferativeTypesCountWriter.hpp"
#include "CellProliferativeTypesWriter.hpp"
#include "CellMutationStatesWriter.hpp"
#include "CellVolumesWriter.hpp"
#include "CellIdWriter.hpp"
#include "Warnings.hpp"
#include "AbstractCellBasedTestSuite.hpp"
#include "SmartPointers.hpp"
#include "PetscSetupAndFinalize.hpp"
#include "Debug.hpp"

Running Multiple Simulations# First of all, we define the test class.

class TestPottsCrypt : public AbstractCellBasedTestSuite
{
public :
void TestMultiplePottsCrypt () throw ( Exception )
{

These numbers specify the random seeds to use

unsigned start_sim = 1 ;
unsigned num_sims = 1 ;

Times to sample the cell number

double mid_time = 25 ;
double end_time = 50 ;

Specifies the alues of

and

to use

double cuberoot10 = 2.15443469 ;
double temp [ 10 ] = { 0.001 , 0.001 * cuberoot10 , 0.001 * cuberoot10 * cuberoot10 , 0.01 , 0.01 * cuberoot10 , 0.01 * cuberoot10 * cuberoot10 , 0.1 , 0.1 * cuberoot10 , 0.1 * cuberoot10 * cuberoot10 , 1.0 }; // 3.1623, 10.0}; // 31.623, 100};
unsigned max_temp_index = 10 ;
double dt [ 7 ] = { 0.1 , 0.01 * cuberoot10 * cuberoot10 , 0.01 * cuberoot10 , 0.01 , 0.001 * cuberoot10 * cuberoot10 , 0.001 * cuberoot10 , 0.001 };
unsigned max_dt_index = 7 ;

This code allows us to to output the number of cells in the crypt, at the mid and end point, for each simulation in the crypt to a file,

This is what’s used to make the contour plots in Figure 4 .

out_stream p_cell_number_file ;
OutputFileHandler output_file_handler ( "Potts/CylindricalCrypt/Sweeps/" , false );
p_cell_number_file = output_file_handler . OpenOutputFile ( "cellnumbers.dat" );

First loop over

.

for ( unsigned dt_index = 0 ; dt_index < max_dt_index ; dt_index ++ )
{
std :: cout << " \n Dt " << dt [ dt_index ] << "... " << std :: flush ;

Then loop over

.

for ( unsigned temp_index = 0 ; temp_index < max_temp_index ; temp_index ++ )
{
std :: cout << " \n\t Temp " << temp [ temp_index ] << ", " << std :: flush ;
double number_of_cells_in_middle = 0.0 ;
double number_of_cells_at_end = 0.0 ;

Finally loop over the random seed.

for ( unsigned index = start_sim ; index < start_sim + num_sims ; index ++ )
{
std :: cout << " Run number " << index << "... " << std :: flush ;

Re seed the random number generator

RandomNumberGenerator :: Instance () -> Reseed ( 100 * index );

THe rest is very simular to the single simulation

double crypt_length = 100 ;

Create a simple 2D PottsMesh and some cells

PottsMeshGenerator < 2 > generator ( 50 , 10 , 5 , 110 , 20 , 5 , 1 , 1 , 1 , true , true );
PottsMesh < 2 >* p_mesh = generator . GetMesh ();
// Create cells
MAKE_PTR ( TransitCellProliferativeType , p_transit_type );
std :: vector < CellPtr > cells ;
CellsGenerator < SimpleWntCellCycleModel , 2 > cells_generator ;
cells_generator . GenerateBasicRandom ( cells , p_mesh -> GetNumElements (), p_transit_type );
// Alter cells properties
for ( unsigned i = 0 ; i < cells . size (); i ++ )
{
dynamic_cast < SimpleWntCellCycleModel *> ( cells [ i ] -> GetCellCycleModel ()) -> SetTransitCellG1Duration ( 6 );
dynamic_cast < SimpleWntCellCycleModel *> ( cells [ i ] -> GetCellCycleModel ()) -> SetWntTransitThreshold ( 2.0 / 3.0 );
}

Create cell population, Wnt stimulus and the simulaton.

PottsBasedCellPopulation < 2 > cell_population ( * p_mesh , cells );
cell_population . SetNumSweepsPerTimestep ( 1 );

Select the appropriate

.

cell_population . SetTemperature ( temp [ temp_index ]);
cell_population . AddCellPopulationCountWriter < CellProliferativeTypesCountWriter > ();
cell_population . AddCellWriter < CellProliferativeTypesWriter > ();
cell_population . AddCellWriter < CellMutationStatesWriter > ();
cell_population . AddCellWriter < CellVolumesWriter > ();
cell_population . AddCellWriter < CellIdWriter > ();
// Create an instance of a Wnt concentration
WntConcentration < 2 >:: Instance () -> SetType ( LINEAR );
WntConcentration < 2 >:: Instance () -> SetCellPopulation ( cell_population );
WntConcentration < 2 >:: Instance () -> SetCryptLength ( crypt_length );
// Set up cell-based simulation
OnLatticeSimulation < 2 > simulator ( cell_population );

Select the appropriate

simulator . SetDt ( dt [ dt_index ]);
simulator . SetSamplingTimestepMultiple (( unsigned )( 1.0 / dt [ dt_index ]));
simulator . SetOutputCellVelocities ( true );

Create output directory, this is based on the loops.

std :: stringstream out ;
out << "/Dt_" << dt [ dt_index ] << "/Temp_" << temp [ temp_index ] << "/RunIndex_" << index ;
std :: string output_directory = "Potts/CylindricalCrypt/Sweeps/" + out . str ();
simulator . SetOutputDirectory ( output_directory );

Create cell killers and update rules.

MAKE_PTR_ARGS ( SloughingCellKiller < 2 > , p_sloughing_killer , ( & cell_population , crypt_length ));
simulator . AddCellKiller ( p_sloughing_killer );
MAKE_PTR ( VolumeConstraintPottsUpdateRule < 2 > , p_volume_constraint_update_rule );
p_volume_constraint_update_rule -> SetMatureCellTargetVolume ( 25 );
p_volume_constraint_update_rule -> SetDeformationEnergyParameter ( 0.1 ); //Default is 0.5
simulator . AddPottsUpdateRule ( p_volume_constraint_update_rule );
MAKE_PTR ( AdhesionPottsUpdateRule < 2 > , p_adhesion_update_rule );
simulator . AddPottsUpdateRule ( p_adhesion_update_rule );

Run simulation to middle, and store the number of cells.

simulator . SetEndTime ( mid_time );
simulator . Solve ();
// Get number of elements of non zero size
PottsMesh < 2 >& potts_mesh = static_cast < PottsMesh < 2 >&> ( simulator . rGetCellPopulation (). rGetMesh ());
unsigned local_num_cells_in_middle = 0 ;
for ( PottsMesh < 2 >:: PottsElementIterator elem_iter = potts_mesh . GetElementIteratorBegin ();
elem_iter != potts_mesh . GetElementIteratorEnd ();
++ elem_iter )
{
if ( elem_iter -> GetNumNodes () > 0 )
{
local_num_cells_in_middle ++ ;
}
}
number_of_cells_in_middle += local_num_cells_in_middle ;

Run simulation to end, and store the number of cells.

simulator . SetEndTime ( end_time );
simulator . Solve ();
// Get number of elements of non zero size
PottsMesh < 2 >& potts_mesh_end = static_cast < PottsMesh < 2 >&> ( simulator . rGetCellPopulation (). rGetMesh ());
unsigned local_num_cells_at_end = 0 ;
for ( PottsMesh < 2 >:: PottsElementIterator elem_iter = potts_mesh_end . GetElementIteratorBegin ();
elem_iter != potts_mesh_end . GetElementIteratorEnd ();
++ elem_iter )
{
if ( elem_iter -> GetNumNodes () > 0 )
{
local_num_cells_at_end ++ ;
}
}
number_of_cells_at_end += local_num_cells_at_end ;

Finally we reset singletons as we’re running multiple simulations in a loop.

WntConcentration < 2 >:: Destroy ();
SimulationTime :: Destroy ();
SimulationTime :: Instance () -> SetStartTime ( 0.0 );
RandomNumberGenerator :: Destroy ();
* p_cell_number_file << temp [ temp_index ] << " \t " << dt [ dt_index ] << " \t " << index << " \t " << local_num_cells_in_middle << " \t " << local_num_cells_at_end << " \n " ;
}
number_of_cells_in_middle = 0.0 ;
number_of_cells_at_end = 0.0 ;
}
}
std :: cout << " \n " << std :: flush ;
p_cell_number_file -> close ();
}

With the parameters as above the simulation will take a couple of hours, this is due to sweeping over very small

s.

The data to reproduce Figure 4 can be generated by running this simulation for more random seeds and averaging the results as described in the paper.
The data is in the

`/tmp/$USER/testoutput/Potts/PottsCryptSweeps/cellnumbers.dat`

file.

Code# The full code is given below

File name `TestPottsCryptSweepsLiteratePaper.hpp`

# #include <cxxtest/TestSuite.h>
#include "CellBasedSimulationArchiver.hpp"
#include "TransitCellProliferativeType.hpp"
#include "PottsMeshGenerator.hpp"
#include "CellsGenerator.hpp"
#include "WntConcentration.hpp"
#include "SimpleWntCellCycleModel.hpp"
#include "PottsBasedCellPopulation.hpp"
#include "VolumeConstraintPottsUpdateRule.hpp"
#include "AdhesionPottsUpdateRule.hpp"
#include "SloughingCellKiller.hpp"
#include "OnLatticeSimulation.hpp"
#include "CellProliferativeTypesCountWriter.hpp"
#include "CellProliferativeTypesWriter.hpp"
#include "CellMutationStatesWriter.hpp"
#include "CellVolumesWriter.hpp"
#include "CellIdWriter.hpp"
#include "Warnings.hpp"
#include "AbstractCellBasedTestSuite.hpp"
#include "SmartPointers.hpp"
#include "PetscSetupAndFinalize.hpp"
#include "Debug.hpp"
class TestPottsCrypt : public AbstractCellBasedTestSuite
{
public :
void TestMultiplePottsCrypt () throw ( Exception )
{
unsigned start_sim = 1 ;
unsigned num_sims = 1 ;
double mid_time = 25 ;
double end_time = 50 ;
double cuberoot10 = 2.15443469 ;
double temp [ 10 ] = { 0.001 , 0.001 * cuberoot10 , 0.001 * cuberoot10 * cuberoot10 , 0.01 , 0.01 * cuberoot10 , 0.01 * cuberoot10 * cuberoot10 , 0.1 , 0.1 * cuberoot10 , 0.1 * cuberoot10 * cuberoot10 , 1.0 }; // 3.1623, 10.0}; // 31.623, 100};
unsigned max_temp_index = 10 ;
double dt [ 7 ] = { 0.1 , 0.01 * cuberoot10 * cuberoot10 , 0.01 * cuberoot10 , 0.01 , 0.001 * cuberoot10 * cuberoot10 , 0.001 * cuberoot10 , 0.001 };
unsigned max_dt_index = 7 ;
out_stream p_cell_number_file ;
OutputFileHandler output_file_handler ( "Potts/CylindricalCrypt/Sweeps/" , false );
p_cell_number_file = output_file_handler . OpenOutputFile ( "cellnumbers.dat" );
for ( unsigned dt_index = 0 ; dt_index < max_dt_index ; dt_index ++ )
{
std :: cout << " \n Dt " << dt [ dt_index ] << "... " << std :: flush ;
for ( unsigned temp_index = 0 ; temp_index < max_temp_index ; temp_index ++ )
{
std :: cout << " \n\t Temp " << temp [ temp_index ] << ", " << std :: flush ;
double number_of_cells_in_middle = 0.0 ;
double number_of_cells_at_end = 0.0 ;
for ( unsigned index = start_sim ; index < start_sim + num_sims ; index ++ )
{
std :: cout << " Run number " << index << "... " << std :: flush ;
RandomNumberGenerator :: Instance () -> Reseed ( 100 * index );
double crypt_length = 100 ;
PottsMeshGenerator < 2 > generator ( 50 , 10 , 5 , 110 , 20 , 5 , 1 , 1 , 1 , true , true );
PottsMesh < 2 >* p_mesh = generator . GetMesh ();
// Create cells
MAKE_PTR ( TransitCellProliferativeType , p_transit_type );
std :: vector < CellPtr > cells ;
CellsGenerator < SimpleWntCellCycleModel , 2 > cells_generator ;
cells_generator . GenerateBasicRandom ( cells , p_mesh -> GetNumElements (), p_transit_type );
// Alter cells properties
for ( unsigned i = 0 ; i < cells . size (); i ++ )
{
dynamic_cast < SimpleWntCellCycleModel *> ( cells [ i ] -> GetCellCycleModel ()) -> SetTransitCellG1Duration ( 6 );
dynamic_cast < SimpleWntCellCycleModel *> ( cells [ i ] -> GetCellCycleModel ()) -> SetWntTransitThreshold ( 2.0 / 3.0 );
}
PottsBasedCellPopulation < 2 > cell_population ( * p_mesh , cells );
cell_population . SetNumSweepsPerTimestep ( 1 );
cell_population . SetTemperature ( temp [ temp_index ]);
cell_population . AddCellPopulationCountWriter < CellProliferativeTypesCountWriter > ();
cell_population . AddCellWriter < CellProliferativeTypesWriter > ();
cell_population . AddCellWriter < CellMutationStatesWriter > ();
cell_population . AddCellWriter < CellVolumesWriter > ();
cell_population . AddCellWriter < CellIdWriter > ();
// Create an instance of a Wnt concentration
WntConcentration < 2 >:: Instance () -> SetType ( LINEAR );
WntConcentration < 2 >:: Instance () -> SetCellPopulation ( cell_population );
WntConcentration < 2 >:: Instance () -> SetCryptLength ( crypt_length );
// Set up cell-based simulation
OnLatticeSimulation < 2 > simulator ( cell_population );
simulator . SetDt ( dt [ dt_index ]);
simulator . SetSamplingTimestepMultiple (( unsigned )( 1.0 / dt [ dt_index ]));
simulator . SetOutputCellVelocities ( true );
std :: stringstream out ;
out << "/Dt_" << dt [ dt_index ] << "/Temp_" << temp [ temp_index ] << "/RunIndex_" << index ;
std :: string output_directory = "Potts/CylindricalCrypt/Sweeps/" + out . str ();
simulator . SetOutputDirectory ( output_directory );
MAKE_PTR_ARGS ( SloughingCellKiller < 2 > , p_sloughing_killer , ( & cell_population , crypt_length ));
simulator . AddCellKiller ( p_sloughing_killer );
MAKE_PTR ( VolumeConstraintPottsUpdateRule < 2 > , p_volume_constraint_update_rule );
p_volume_constraint_update_rule -> SetMatureCellTargetVolume ( 25 );
p_volume_constraint_update_rule -> SetDeformationEnergyParameter ( 0.1 ); //Default is 0.5
simulator . AddPottsUpdateRule ( p_volume_constraint_update_rule );
MAKE_PTR ( AdhesionPottsUpdateRule < 2 > , p_adhesion_update_rule );
simulator . AddPottsUpdateRule ( p_adhesion_update_rule );
simulator . SetEndTime ( mid_time );
simulator . Solve ();
// Get number of elements of non zero size
PottsMesh < 2 >& potts_mesh = static_cast < PottsMesh < 2 >&> ( simulator . rGetCellPopulation (). rGetMesh ());
unsigned local_num_cells_in_middle = 0 ;
for ( PottsMesh < 2 >:: PottsElementIterator elem_iter = potts_mesh . GetElementIteratorBegin ();
elem_iter != potts_mesh . GetElementIteratorEnd ();
++ elem_iter )
{
if ( elem_iter -> GetNumNodes () > 0 )
{
local_num_cells_in_middle ++ ;
}
}
number_of_cells_in_middle += local_num_cells_in_middle ;
simulator . SetEndTime ( end_time );
simulator . Solve ();
// Get number of elements of non zero size
PottsMesh < 2 >& potts_mesh_end = static_cast < PottsMesh < 2 >&> ( simulator . rGetCellPopulation (). rGetMesh ());
unsigned local_num_cells_at_end = 0 ;
for ( PottsMesh < 2 >:: PottsElementIterator elem_iter = potts_mesh_end . GetElementIteratorBegin ();
elem_iter != potts_mesh_end . GetElementIteratorEnd ();
++ elem_iter )
{
if ( elem_iter -> GetNumNodes () > 0 )
{
local_num_cells_at_end ++ ;
}
}
number_of_cells_at_end += local_num_cells_at_end ;
WntConcentration < 2 >:: Destroy ();
SimulationTime :: Destroy ();
SimulationTime :: Instance () -> SetStartTime ( 0.0 );
RandomNumberGenerator :: Destroy ();
* p_cell_number_file << temp [ temp_index ] << " \t " << dt [ dt_index ] << " \t " << index << " \t " << local_num_cells_in_middle << " \t " << local_num_cells_at_end << " \n " ;
}
number_of_cells_in_middle = 0.0 ;
number_of_cells_at_end = 0.0 ;
}
}
std :: cout << " \n " << std :: flush ;
p_cell_number_file -> close ();
}
};