This tutorial was generated from the file projects/CryptProliferationDistribution/test/TestRecoveredCryptLiteratePaper.hpp at revision r25342. 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 "CellProliferativeTypesCountWriter.hpp"
#include "CellMutationStatesCountWriter.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 TestRecoveredCryptLiteratePaper : 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)
    {

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

//    	double percent_mutant = 0.3;

		bool contact_inhibition = true;
		double end_time = 5200;

    	// Optimal Healthy Model
    	unsigned healthy_cell_proliferation_model = 3u; // Spatial Wnt at birth
    	bool healthy_wnt_dependend_ccd = true;
		double healthy_wnt_thresh = 0.6;
		double healthy_CI = 0.9;

		//Optimal Mutant model, as above with
		unsigned mutant_cell_proliferation_model = 3u; // Spatial Wnt at birth
		bool mutant_wnt_dependend_ccd = true;
		assert(mutant_wnt_dependend_ccd == healthy_wnt_dependend_ccd); // Currently needs to be the same in the CCM
	    double mutant_wnt_thresh = 0.5;
		double mutant_CI = 0.6;

        // 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_percents;   index ++)
//         {
//        	double percent_mutant = min_percent + (double)index / double(num_percents) * (max_percent-min_percent);
//
			PRINT_VARIABLE(percent_mutant);

			// 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);
			boost::shared_ptr<AbstractCellProperty> p_mutant_state(CellPropertyRegistry::Instance()->Get<ApcTwoHitCellMutationState>());
			boost::shared_ptr<AbstractCellProperty> p_paneth_mutant_state(CellPropertyRegistry::Instance()->Get<ApcOneHitCellMutationState>());

			// 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(healthy_cell_proliferation_model);
				  dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetIsContactInhibitionCellCycleDuration((bool)contact_inhibition);
				  dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetIsWntDependentCellCycleDuration((bool)healthy_wnt_dependend_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 (healthy_cell_proliferation_model ==1 ) // i.e Pedigree dependent
				  {
					  assert(0);
//				    	  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(healthy_wnt_thresh);
					  dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetMutantWntThreshold(mutant_wnt_thresh);
					  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(healthy_CI);
				  dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetMutantQuiescentVolumeFraction(mutant_CI);

mutant cells

				  dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetMutantHealthyRatio(percent_mutant);
				  dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetHealthyCellProliferationModel(healthy_cell_proliferation_model);
				  dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetMutantCellProliferationModel(mutant_cell_proliferation_model);

			}

			// 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);
					// Give them a mutation to make them easier to track (doesn't actually do anything but label the cells)
				    cells[index]->SetMutationState(p_paneth_mutant_state);
			}

			// Make some cells mutant
			for(unsigned cell_index= 0;  cell_index<num_cells; cell_index++)
			{

Generate a uniform random number to choose between Healthy and mutant cell in appropriate ratio Contact Inhibition specific parameters

				if(!cells[cell_index]->GetCellProliferativeType()->IsType<PanethCellProliferativeType>())
				{
					double u = RandomNumberGenerator::Instance()->ranf();
					//PRINT_2_VARIABLES(u,percent_mutant);
					if (u < percent_mutant)// Mutant cell
					{
						cells[cell_index]->SetMutationState(p_mutant_state);
						//dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetCellProliferationModel(4u);
					}
				}
			}

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

			// Output data
			crypt.AddCellPopulationCountWriter<CellProliferativeTypesCountWriter>();
			crypt.AddCellPopulationCountWriter<CellMutationStatesCountWriter>();
			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_volume_modifier);
			simulator.AddSimulationModifier(p_volume_modifier);

			//Create output directory
			std::stringstream out;
			if (healthy_cell_proliferation_model == 1)
			{
				assert(0);
				//out << "FitCCM_"<< cell_proliferation_model << "_CI_" << contact_inhibition << "_WDCCD_" << wnt_dependent_ccd << "_MaxGen_" << param << "_CIthresh_" << CIparam;
			}
			else
			{
				out << "Fit_PercentMutant_" << percent_mutant;
			}
			std::string output_directory = "CryptMutantFit/" +  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 TestRecoveredCryptLiteratePaper.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 "CellProliferativeTypesCountWriter.hpp"
#include "CellMutationStatesCountWriter.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 TestRecoveredCryptLiteratePaper : 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)
    {

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

//    	double percent_mutant = 0.3;

		bool contact_inhibition = true;
		double end_time = 5200;

    	// Optimal Healthy Model
    	unsigned healthy_cell_proliferation_model = 3u; // Spatial Wnt at birth
    	bool healthy_wnt_dependend_ccd = true;
		double healthy_wnt_thresh = 0.6;
		double healthy_CI = 0.9;

		//Optimal Mutant model, as above with
		unsigned mutant_cell_proliferation_model = 3u; // Spatial Wnt at birth
		bool mutant_wnt_dependend_ccd = true;
		assert(mutant_wnt_dependend_ccd == healthy_wnt_dependend_ccd); // Currently needs to be the same in the CCM
	    double mutant_wnt_thresh = 0.5;
		double mutant_CI = 0.6;

        // 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_percents;   index ++)
//         {
//        	double percent_mutant = min_percent + (double)index / double(num_percents) * (max_percent-min_percent);
//
			PRINT_VARIABLE(percent_mutant);

			// 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);
			boost::shared_ptr<AbstractCellProperty> p_mutant_state(CellPropertyRegistry::Instance()->Get<ApcTwoHitCellMutationState>());
			boost::shared_ptr<AbstractCellProperty> p_paneth_mutant_state(CellPropertyRegistry::Instance()->Get<ApcOneHitCellMutationState>());

			// 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);
				  dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetCellProliferationModel(healthy_cell_proliferation_model);
				  dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetIsContactInhibitionCellCycleDuration((bool)contact_inhibition);
				  dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetIsWntDependentCellCycleDuration((bool)healthy_wnt_dependend_ccd);

				  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

				  if (healthy_cell_proliferation_model ==1 ) // i.e Pedigree dependent
				  {
					  assert(0);
//				    	  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(healthy_wnt_thresh);
					  dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetMutantWntThreshold(mutant_wnt_thresh);
					  dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetMaxTransitGenerations(UINT_MAX);
				  }

				  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(healthy_CI);
				  dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetMutantQuiescentVolumeFraction(mutant_CI);

				  dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetMutantHealthyRatio(percent_mutant);
				  dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetHealthyCellProliferationModel(healthy_cell_proliferation_model);
				  dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetMutantCellProliferationModel(mutant_cell_proliferation_model);

			}

			// 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);
					// Give them a mutation to make them easier to track (doesn't actually do anything but label the cells)
				    cells[index]->SetMutationState(p_paneth_mutant_state);
			}

			// Make some cells mutant
			for(unsigned cell_index= 0;  cell_index<num_cells; cell_index++)
			{
				if(!cells[cell_index]->GetCellProliferativeType()->IsType<PanethCellProliferativeType>())
				{
					double u = RandomNumberGenerator::Instance()->ranf();
					//PRINT_2_VARIABLES(u,percent_mutant);
					if (u < percent_mutant)// Mutant cell
					{
						cells[cell_index]->SetMutationState(p_mutant_state);
						//dynamic_cast<CryptCellCycleModel*>(cells[cell_index]->GetCellCycleModel())->SetCellProliferationModel(4u);
					}
				}
			}

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

			// Output data
			crypt.AddCellPopulationCountWriter<CellProliferativeTypesCountWriter>();
			crypt.AddCellPopulationCountWriter<CellMutationStatesCountWriter>();
			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_volume_modifier);
			simulator.AddSimulationModifier(p_volume_modifier);

			//Create output directory
			std::stringstream out;
			if (healthy_cell_proliferation_model == 1)
			{
				assert(0);
				//out << "FitCCM_"<< cell_proliferation_model << "_CI_" << contact_inhibition << "_WDCCD_" << wnt_dependent_ccd << "_MaxGen_" << param << "_CIthresh_" << CIparam;
			}
			else
			{
				out << "Fit_PercentMutant_" << percent_mutant;
			}
			std::string output_directory = "CryptMutantFit/" +  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);
		}

//	}
};
Prev
Next