| Chaste

This tutorial was generated from the file projects/CryptProliferationDistribution/test/TestCryptProliferationDistributionLiteratePaper.hpp at revision r25435. Note that the code is given in full at the bottom of the page.

#include <cxxtest/TestSuite.h>

// Must be included before any other cell_based headers
#include "CellBasedSimulationArchiver.hpp"
#include "AbstractCellBasedTestSuite.hpp"
#include "CellBasedEventHandler.hpp"

#include "OffLatticeSimulation.hpp"
#include "VolumeTrackingModifier.hpp"

#include "CryptCellsGenerator.hpp"
#include "CellsGenerator.hpp"

#include "CylindricalHoneycombMeshGenerator.hpp"

#include "CryptCellCycleModel.hpp"

#include "PanethCellProliferativeType.hpp"
#include "TransitCellProliferativeType.hpp"

#include "CellAgesWriter.hpp"
#include "CellVolumesWriter.hpp"
#include "CellProliferativeTypesWriter.hpp"
#include "CellMutationStatesWriter.hpp"
#include "NodeVelocityWriter.hpp"

#include "GeneralisedLinearSpringForce.hpp"
#include "DifferentialAdhesionSpringForce.hpp"
#include "RepulsionForce.hpp"
#include "CellRetainerForce.hpp"

#include "CryptSimulationBoundaryCondition.hpp"
#include "CryptGeometryBoundaryCondition3d.hpp"

#include "MeshBasedCellPopulationWithGhostNodes.hpp"
#include "NodeBasedCellPopulation.hpp"

#include "SloughingCellKiller.hpp"
#include "PlaneBasedCellKiller.hpp"

#include "PetscSetupAndFinalize.hpp"

#include "Debug.hpp"

class TestCryptProliferationDistributionLiteratePaper : public AbstractCellBasedTestSuite
{
private:

    double mLastStartTime;
    void setUp()
    {
        mLastStartTime = std::clock();
        AbstractCellBasedTestSuite::setUp();
    }
    void tearDown()
    {
        double time = std::clock();
        double elapsed_time = (time - mLastStartTime)/(CLOCKS_PER_SEC);
        std::cout << "Elapsed time: " << elapsed_time << std::endl;
        AbstractCellBasedTestSuite::tearDown();
    }

public:

    void Test3DCrypt() throw (Exception)
    {
    	// Model setup
        // 1 - Pedigree
        // 2 - Spatial Wnt
        // 3 - Spatial Wnt at birth
        // 4 - Mutant

    	TS_ASSERT(CommandLineArguments::Instance()->OptionExists("-CCM"));
    	unsigned cell_proliferation_model = (unsigned) atof(CommandLineArguments::Instance()->GetStringCorrespondingToOption("-CCM").c_str());
    	assert(cell_proliferation_model==1 || cell_proliferation_model==2 || cell_proliferation_model==3  || cell_proliferation_model==4);

    	bool contact_inhibition = CommandLineArguments::Instance()->OptionExists("-CI");

    	bool wnt_dependent_ccd = CommandLineArguments::Instance()->OptionExists("-WDCCD");

    	// Sim and sweep Params
    	TS_ASSERT(CommandLineArguments::Instance()->OptionExists("-end_time"));
    	double end_time = atof(CommandLineArguments::Instance()->GetStringCorrespondingToOption("-end_time").c_str());

    	TS_ASSERT(CommandLineArguments::Instance()->OptionExists("-min"));
        double min_param = atof(CommandLineArguments::Instance()->GetStringCorrespondingToOption("-min").c_str());

        TS_ASSERT(CommandLineArguments::Instance()->OptionExists("-max"));
        double max_param = atof(CommandLineArguments::Instance()->GetStringCorrespondingToOption("-max").c_str());

        assert(min_param>=0);
        if (cell_proliferation_model ==2 || cell_proliferation_model==3) // i.e wnt dependent
        {
        	assert (max_param<=1);
		}

        TS_ASSERT(CommandLineArguments::Instance()->OptionExists("-num_sweeps"));
    	double num_sweeps = atof(CommandLineArguments::Instance()->GetStringCorrespondingToOption("-num_sweeps").c_str());

    	TS_ASSERT(CommandLineArguments::Instance()->OptionExists("-min_CI"));
        double min_CI_param = atof(CommandLineArguments::Instance()->GetStringCorrespondingToOption("-min_CI").c_str());

        TS_ASSERT(CommandLineArguments::Instance()->OptionExists("-max_CI"));
        double max_CI_param = atof(CommandLineArguments::Instance()->GetStringCorrespondingToOption("-max_CI").c_str());

    	TS_ASSERT(CommandLineArguments::Instance()->OptionExists("-num_CI_sweeps"));
    	double num_CI_sweeps = atof(CommandLineArguments::Instance()->GetStringCorrespondingToOption("-num_CI_sweeps").c_str());

        // Crypt Setup
        double cell_radius = 3.5;
        double crypt_length = 70; //70
        double crypt_radius = 8.0/M_PI*6.0;                 // Choosing same dimensions as for halted migration paper
        // For this size domain there are about 75 cells in the bottom hemisphere so this makes about 20% Paneth cells

		unsigned num_paneth_cells = 15;//15;
		unsigned num_stem_cells = 60;
		unsigned num_cells = num_paneth_cells + num_stem_cells;

        double stem_retainer_force_magnitude = 7.5*10;
        double paneth_retainer_force_magnitude = 7.5*10;

        for (unsigned index = 0; index<=num_sweeps;   index ++)
        {
        	double param = min_param + (double)index / double(num_sweeps) * (max_param-min_param);

			PRINT_3_VARIABLES(cell_proliferation_model,
								index,
								param);

			// Extra code to calculate the generation if Model 1
			unsigned upper_max_transit_generations = (unsigned)ceil(param);
			unsigned lower_max_transit_generations = (unsigned)floor(param);
			double prob_of_upper_max_transit_generations = param - floor(param);

			//Sweep over CI param
			for (unsigned CIindex = 0; CIindex<=num_CI_sweeps;   CIindex ++)
			{
				double CIparam = min_CI_param + (double)CIindex / double(num_CI_sweeps) * (max_CI_param-min_CI_param);
				PRINT_2_VARIABLES(CIindex,CIparam);

				// Create some starter nodes
				std::vector<Node<3>*> nodes;
				for(unsigned node_index= 0;  node_index<num_cells; node_index++)
				{
					double x = crypt_radius/2.0 * sin(node_index*2.0*M_PI/num_cells);
					double y = crypt_radius/2.0 * cos(node_index*2.0*M_PI/num_cells);
					double z = 0.0;
					nodes.push_back(new Node<3>(node_index, false, x, y, z));
				}

				// Convert this to a NodesOnlyMesh
				NodesOnlyMesh<3> mesh;
				mesh.ConstructNodesWithoutMesh(nodes,cell_radius*3.0);

				MAKE_PTR(PanethCellProliferativeType, p_paneth_type);

				// Create cells
				std::vector<CellPtr> cells;

				CellsGenerator<CryptCellCycleModel, 3> cells_generator;
				cells_generator.GenerateBasicRandom(cells, mesh.GetNumNodes());

				//Change properties of the ccm
				for (unsigned cell_index= 0;  cell_index<cells.size(); cell_index++)
				{
					cells[cell_index]->GetCellData()->SetItem("Radius", cell_radius);

					// Specify CCM
					dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetCellProliferationModel(cell_proliferation_model);
					dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetIsContactInhibitionCellCycleDuration((bool)contact_inhibition);
					dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetIsWntDependentCellCycleDuration((bool)wnt_dependent_ccd);

                    // Set some default CCD parameters So total CCM is U[10,14] and (U[22,26] at base if variable)
                    dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetMDuration(4.0);
                    dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetSDuration(4.0);
                    dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetG2Duration(2.0);
                    dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetTransitCellG1Duration(2.0);  // so total CCM is U[10,14] at threshold
                    dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetStemCellG1Duration(14.0);  // so total CCM is U[10,14] at base

					//Threshold and Generation specific parameters
					if (cell_proliferation_model ==1 ) // i.e Pedigree dependent
				    {
				    	  dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetWntThreshold(1.0);
				    	  if (RandomNumberGenerator::Instance()->ranf()<prob_of_upper_max_transit_generations)
				    	  {
				    		  dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetMaxTransitGenerations(upper_max_transit_generations); // Mutant = MAX_UNSIGNED
				    	  }
				    	  else
				    	  {
				    		  dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetMaxTransitGenerations(lower_max_transit_generations); // Mutant = MAX_UNSIGNED
				    	  }
				    }
				    else
				    {
				    	  dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetWntThreshold(param);
				    	  dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetMaxTransitGenerations(UINT_MAX);
  				    }

                    // Contact Inhibition specific parameters (Mutant, same CI)
                    dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetEquilibriumVolume(M_PI*4.0/3.0*cell_radius*cell_radius*cell_radius);
                    dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetQuiescentVolumeFraction(CIparam);
                    dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetMutantQuiescentVolumeFraction(CIparam);

				}

				// Make some cells paneth cells
				for(unsigned cell_index= 0;  cell_index<num_paneth_cells; cell_index++)
				{
                    unsigned index = cell_index * num_cells/num_paneth_cells;
                    if (index > num_cells)
                    {
                        index = num_cells;
                    }
                    cells[index]->SetCellProliferativeType(p_paneth_type);
				}

				// Create cell population
				NodeBasedCellPopulation<3> crypt(mesh, cells);
				crypt.SetUseVariableRadii(true);

				// Output data
				crypt.AddCellWriter<CellAgesWriter>();
				crypt.AddCellWriter<CellVolumesWriter>();
				crypt.AddCellWriter<CellProliferativeTypesWriter>();
				crypt.AddCellWriter<CellMutationStatesWriter>();
				crypt.AddPopulationWriter<NodeVelocityWriter>();

				crypt.SetAbsoluteMovementThreshold(50.0);

				// Create an instance of a Wnt concentration NOTE DO THIS BEFORE THE SIMULATION OTHERWISE CELLS CANT INITIALISE
				WntConcentration<3>::Instance()->SetType(LINEAR);
				WntConcentration<3>::Instance()->SetCellPopulation(crypt);
				WntConcentration<3>::Instance()->SetCryptLength(crypt_length);

				// Create a contact inhibition simulator
				OffLatticeSimulation<3> simulator(crypt);

				simulator.SetOutputDivisionLocations(true);
				simulator.SetDt(1.0/200.0);
				simulator.SetSamplingTimestepMultiple(200);
				simulator.SetEndTime(end_time);

				// Add Volume Tracking Modifier
				MAKE_PTR(VolumeTrackingModifier<3>, p_modifier);
				simulator.AddSimulationModifier(p_modifier);

				//Create output directory
				std::stringstream out;
				if (cell_proliferation_model == 1)
				{
					out << "FitCCM_"<< cell_proliferation_model << "_CI_" << contact_inhibition << "_WDCCD_" << wnt_dependent_ccd << "_MaxGen_" << param << "_CIthresh_" << CIparam;

				}
				else
				{
					out << "FitCCM_"<< cell_proliferation_model << "_CI_" << contact_inhibition << "_WDCCD_" << wnt_dependent_ccd << "_WntThresh_" << param << "_CIthresh_" << CIparam;
				}
				std::string output_directory = "CryptProlifFit/" +  out.str();
				simulator.SetOutputDirectory(output_directory);

				// Create a force law and pass it to the simulation
				MAKE_PTR(DifferentialAdhesionSpringForce<3>, p_force);
				p_force->SetMeinekeSpringStiffness(30.0); //normally 15.0 but 30 in all CellBased Papers;
				p_force->SetCutOffLength(cell_radius*3.0);
				simulator.AddForce(p_force);

				// Apply a retainer to keep stem and paneth cells at the base of the crypt
				MAKE_PTR(CellRetainerForce<3>, p_retainer_force);
				p_retainer_force->SetStemCellForceMagnitudeParameter(stem_retainer_force_magnitude);
				p_retainer_force->SetPanethCellForceMagnitudeParameter(paneth_retainer_force_magnitude);
				simulator.AddForce(p_retainer_force);

				// Apply a boundary condition to represent a 3d crypt
				MAKE_PTR_ARGS(CryptGeometryBoundaryCondition3d<3>, p_boundary_condition, (&crypt, 0.0));
				simulator.AddCellPopulationBoundaryCondition(p_boundary_condition);

				// Create cell killer and pass in to crypt simulation
				MAKE_PTR_ARGS(PlaneBasedCellKiller<3>, p_cell_killer,(&crypt, crypt_length*unit_vector<double>(3,2), unit_vector<double>(3,2)));
				simulator.AddCellKiller(p_cell_killer);

				// Run simulation
				simulator.Solve();

				// Extra Gubbins to get to loop: this is usually done by the SetUp and TearDown methods
				WntConcentration<3>::Instance()->Destroy();
				SimulationTime::Instance()->Destroy();
				SimulationTime::Instance()->SetStartTime(0.0);
			}
         }
	}
};

Code

The full code is given below

File name `TestCryptProliferationDistributionLiteratePaper.hpp`

#include <cxxtest/TestSuite.h>

// Must be included before any other cell_based headers
#include "CellBasedSimulationArchiver.hpp"
#include "AbstractCellBasedTestSuite.hpp"
#include "CellBasedEventHandler.hpp"

#include "OffLatticeSimulation.hpp"
#include "VolumeTrackingModifier.hpp"

#include "CryptCellsGenerator.hpp"
#include "CellsGenerator.hpp"

#include "CylindricalHoneycombMeshGenerator.hpp"

#include "CryptCellCycleModel.hpp"

#include "PanethCellProliferativeType.hpp"
#include "TransitCellProliferativeType.hpp"

#include "CellAgesWriter.hpp"
#include "CellVolumesWriter.hpp"
#include "CellProliferativeTypesWriter.hpp"
#include "CellMutationStatesWriter.hpp"
#include "NodeVelocityWriter.hpp"

#include "GeneralisedLinearSpringForce.hpp"
#include "DifferentialAdhesionSpringForce.hpp"
#include "RepulsionForce.hpp"
#include "CellRetainerForce.hpp"

#include "CryptSimulationBoundaryCondition.hpp"
#include "CryptGeometryBoundaryCondition3d.hpp"

#include "MeshBasedCellPopulationWithGhostNodes.hpp"
#include "NodeBasedCellPopulation.hpp"

#include "SloughingCellKiller.hpp"
#include "PlaneBasedCellKiller.hpp"

#include "PetscSetupAndFinalize.hpp"

#include "Debug.hpp"

class TestCryptProliferationDistributionLiteratePaper : public AbstractCellBasedTestSuite
{
private:

    double mLastStartTime;
    void setUp()
    {
        mLastStartTime = std::clock();
        AbstractCellBasedTestSuite::setUp();
    }
    void tearDown()
    {
        double time = std::clock();
        double elapsed_time = (time - mLastStartTime)/(CLOCKS_PER_SEC);
        std::cout << "Elapsed time: " << elapsed_time << std::endl;
        AbstractCellBasedTestSuite::tearDown();
    }

public:

    void Test3DCrypt() throw (Exception)
    {
    	// Model setup
        // 1 - Pedigree
        // 2 - Spatial Wnt
        // 3 - Spatial Wnt at birth
        // 4 - Mutant

    	TS_ASSERT(CommandLineArguments::Instance()->OptionExists("-CCM"));
    	unsigned cell_proliferation_model = (unsigned) atof(CommandLineArguments::Instance()->GetStringCorrespondingToOption("-CCM").c_str());
    	assert(cell_proliferation_model==1 || cell_proliferation_model==2 || cell_proliferation_model==3  || cell_proliferation_model==4);

    	bool contact_inhibition = CommandLineArguments::Instance()->OptionExists("-CI");

    	bool wnt_dependent_ccd = CommandLineArguments::Instance()->OptionExists("-WDCCD");

    	// Sim and sweep Params
    	TS_ASSERT(CommandLineArguments::Instance()->OptionExists("-end_time"));
    	double end_time = atof(CommandLineArguments::Instance()->GetStringCorrespondingToOption("-end_time").c_str());

    	TS_ASSERT(CommandLineArguments::Instance()->OptionExists("-min"));
        double min_param = atof(CommandLineArguments::Instance()->GetStringCorrespondingToOption("-min").c_str());

        TS_ASSERT(CommandLineArguments::Instance()->OptionExists("-max"));
        double max_param = atof(CommandLineArguments::Instance()->GetStringCorrespondingToOption("-max").c_str());

        assert(min_param>=0);
        if (cell_proliferation_model ==2 || cell_proliferation_model==3) // i.e wnt dependent
        {
        	assert (max_param<=1);
		}

        TS_ASSERT(CommandLineArguments::Instance()->OptionExists("-num_sweeps"));
    	double num_sweeps = atof(CommandLineArguments::Instance()->GetStringCorrespondingToOption("-num_sweeps").c_str());

    	TS_ASSERT(CommandLineArguments::Instance()->OptionExists("-min_CI"));
        double min_CI_param = atof(CommandLineArguments::Instance()->GetStringCorrespondingToOption("-min_CI").c_str());

        TS_ASSERT(CommandLineArguments::Instance()->OptionExists("-max_CI"));
        double max_CI_param = atof(CommandLineArguments::Instance()->GetStringCorrespondingToOption("-max_CI").c_str());

    	TS_ASSERT(CommandLineArguments::Instance()->OptionExists("-num_CI_sweeps"));
    	double num_CI_sweeps = atof(CommandLineArguments::Instance()->GetStringCorrespondingToOption("-num_CI_sweeps").c_str());

        // Crypt Setup
        double cell_radius = 3.5;
        double crypt_length = 70; //70
        double crypt_radius = 8.0/M_PI*6.0;                 // Choosing same dimensions as for halted migration paper
        // For this size domain there are about 75 cells in the bottom hemisphere so this makes about 20% Paneth cells

		unsigned num_paneth_cells = 15;//15;
		unsigned num_stem_cells = 60;
		unsigned num_cells = num_paneth_cells + num_stem_cells;

        double stem_retainer_force_magnitude = 7.5*10;
        double paneth_retainer_force_magnitude = 7.5*10;

        for (unsigned index = 0; index<=num_sweeps;   index ++)
        {
        	double param = min_param + (double)index / double(num_sweeps) * (max_param-min_param);

			PRINT_3_VARIABLES(cell_proliferation_model,
								index,
								param);

			// Extra code to calculate the generation if Model 1
			unsigned upper_max_transit_generations = (unsigned)ceil(param);
			unsigned lower_max_transit_generations = (unsigned)floor(param);
			double prob_of_upper_max_transit_generations = param - floor(param);

			//Sweep over CI param
			for (unsigned CIindex = 0; CIindex<=num_CI_sweeps;   CIindex ++)
			{
				double CIparam = min_CI_param + (double)CIindex / double(num_CI_sweeps) * (max_CI_param-min_CI_param);
				PRINT_2_VARIABLES(CIindex,CIparam);

				// Create some starter nodes
				std::vector<Node<3>*> nodes;
				for(unsigned node_index= 0;  node_index<num_cells; node_index++)
				{
					double x = crypt_radius/2.0 * sin(node_index*2.0*M_PI/num_cells);
					double y = crypt_radius/2.0 * cos(node_index*2.0*M_PI/num_cells);
					double z = 0.0;
					nodes.push_back(new Node<3>(node_index, false, x, y, z));
				}

				// Convert this to a NodesOnlyMesh
				NodesOnlyMesh<3> mesh;
				mesh.ConstructNodesWithoutMesh(nodes,cell_radius*3.0);

				MAKE_PTR(PanethCellProliferativeType, p_paneth_type);

				// Create cells
				std::vector<CellPtr> cells;

				CellsGenerator<CryptCellCycleModel, 3> cells_generator;
				cells_generator.GenerateBasicRandom(cells, mesh.GetNumNodes());

				//Change properties of the ccm
				for (unsigned cell_index= 0;  cell_index<cells.size(); cell_index++)
				{
					cells[cell_index]->GetCellData()->SetItem("Radius", cell_radius);

					// Specify CCM
					dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetCellProliferationModel(cell_proliferation_model);
					dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetIsContactInhibitionCellCycleDuration((bool)contact_inhibition);
					dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetIsWntDependentCellCycleDuration((bool)wnt_dependent_ccd);

                    // Set some default CCD parameters So total CCM is U[10,14] and (U[22,26] at base if variable)
                    dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetMDuration(4.0);
                    dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetSDuration(4.0);
                    dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetG2Duration(2.0);
                    dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetTransitCellG1Duration(2.0);  // so total CCM is U[10,14] at threshold
                    dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetStemCellG1Duration(14.0);  // so total CCM is U[10,14] at base

					//Threshold and Generation specific parameters
					if (cell_proliferation_model ==1 ) // i.e Pedigree dependent
				    {
				    	  dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetWntThreshold(1.0);
				    	  if (RandomNumberGenerator::Instance()->ranf()<prob_of_upper_max_transit_generations)
				    	  {
				    		  dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetMaxTransitGenerations(upper_max_transit_generations); // Mutant = MAX_UNSIGNED
				    	  }
				    	  else
				    	  {
				    		  dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetMaxTransitGenerations(lower_max_transit_generations); // Mutant = MAX_UNSIGNED
				    	  }
				    }
				    else
				    {
				    	  dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetWntThreshold(param);
				    	  dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetMaxTransitGenerations(UINT_MAX);
  				    }

                    // Contact Inhibition specific parameters (Mutant, same CI)
                    dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetEquilibriumVolume(M_PI*4.0/3.0*cell_radius*cell_radius*cell_radius);
                    dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetQuiescentVolumeFraction(CIparam);
                    dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetMutantQuiescentVolumeFraction(CIparam);

				}

				// Make some cells paneth cells
				for(unsigned cell_index= 0;  cell_index<num_paneth_cells; cell_index++)
				{
                    unsigned index = cell_index * num_cells/num_paneth_cells;
                    if (index > num_cells)
                    {
                        index = num_cells;
                    }
                    cells[index]->SetCellProliferativeType(p_paneth_type);
				}

				// Create cell population
				NodeBasedCellPopulation<3> crypt(mesh, cells);
				crypt.SetUseVariableRadii(true);

				// Output data
				crypt.AddCellWriter<CellAgesWriter>();
				crypt.AddCellWriter<CellVolumesWriter>();
				crypt.AddCellWriter<CellProliferativeTypesWriter>();
				crypt.AddCellWriter<CellMutationStatesWriter>();
				crypt.AddPopulationWriter<NodeVelocityWriter>();

				crypt.SetAbsoluteMovementThreshold(50.0);

				// Create an instance of a Wnt concentration NOTE DO THIS BEFORE THE SIMULATION OTHERWISE CELLS CANT INITIALISE
				WntConcentration<3>::Instance()->SetType(LINEAR);
				WntConcentration<3>::Instance()->SetCellPopulation(crypt);
				WntConcentration<3>::Instance()->SetCryptLength(crypt_length);

				// Create a contact inhibition simulator
				OffLatticeSimulation<3> simulator(crypt);

				simulator.SetOutputDivisionLocations(true);
				simulator.SetDt(1.0/200.0);
				simulator.SetSamplingTimestepMultiple(200);
				simulator.SetEndTime(end_time);

				// Add Volume Tracking Modifier
				MAKE_PTR(VolumeTrackingModifier<3>, p_modifier);
				simulator.AddSimulationModifier(p_modifier);

				//Create output directory
				std::stringstream out;
				if (cell_proliferation_model == 1)
				{
					out << "FitCCM_"<< cell_proliferation_model << "_CI_" << contact_inhibition << "_WDCCD_" << wnt_dependent_ccd << "_MaxGen_" << param << "_CIthresh_" << CIparam;

				}
				else
				{
					out << "FitCCM_"<< cell_proliferation_model << "_CI_" << contact_inhibition << "_WDCCD_" << wnt_dependent_ccd << "_WntThresh_" << param << "_CIthresh_" << CIparam;
				}
				std::string output_directory = "CryptProlifFit/" +  out.str();
				simulator.SetOutputDirectory(output_directory);

				// Create a force law and pass it to the simulation
				MAKE_PTR(DifferentialAdhesionSpringForce<3>, p_force);
				p_force->SetMeinekeSpringStiffness(30.0); //normally 15.0 but 30 in all CellBased Papers;
				p_force->SetCutOffLength(cell_radius*3.0);
				simulator.AddForce(p_force);

				// Apply a retainer to keep stem and paneth cells at the base of the crypt
				MAKE_PTR(CellRetainerForce<3>, p_retainer_force);
				p_retainer_force->SetStemCellForceMagnitudeParameter(stem_retainer_force_magnitude);
				p_retainer_force->SetPanethCellForceMagnitudeParameter(paneth_retainer_force_magnitude);
				simulator.AddForce(p_retainer_force);

				// Apply a boundary condition to represent a 3d crypt
				MAKE_PTR_ARGS(CryptGeometryBoundaryCondition3d<3>, p_boundary_condition, (&crypt, 0.0));
				simulator.AddCellPopulationBoundaryCondition(p_boundary_condition);

				// Create cell killer and pass in to crypt simulation
				MAKE_PTR_ARGS(PlaneBasedCellKiller<3>, p_cell_killer,(&crypt, crypt_length*unit_vector<double>(3,2), unit_vector<double>(3,2)));
				simulator.AddCellKiller(p_cell_killer);

				// Run simulation
				simulator.Solve();

				// Extra Gubbins to get to loop: this is usually done by the SetUp and TearDown methods
				WntConcentration<3>::Instance()->Destroy();
				SimulationTime::Instance()->Destroy();
				SimulationTime::Instance()->SetStartTime(0.0);
			}
         }
	}
};

Code#

File name TestCryptProliferationDistributionLiteratePaper.hpp#

Code

File name `TestCryptProliferationDistributionLiteratePaper.hpp`