diff --git a/test/Test05CellSorting.hpp b/test/Test05CellSorting.hpp index 7f58035..59c1f6a 100644 --- a/test/Test05CellSorting.hpp +++ b/test/Test05CellSorting.hpp @@ -42,23 +42,412 @@ OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #include "CellBasedSimulationArchiver.hpp" #include "AbstractCellBasedWithTimingsTestSuite.hpp" +#include "CellLabel.hpp" +#include "SmartPointers.hpp" +#include "CellsGenerator.hpp" +#include "UniformG1GenerationalCellCycleModel.hpp" +#include "TransitCellProliferativeType.hpp" +#include "DifferentiatedCellProliferativeType.hpp" + +#include "HeterotypicBoundaryLengthWriter.hpp" + +#include "OffLatticeSimulation.hpp" +#include "VertexBasedCellPopulation.hpp" +#include "HoneycombVertexMeshGenerator.hpp" +#include "NagaiHondaDifferentialAdhesionForce.hpp" +#include "RandomMotionForce.hpp" +#include "SimpleTargetAreaModifier.hpp" + +#include "MeshBasedCellPopulationWithGhostNodes.hpp" +#include "HoneycombMeshGenerator.hpp" +#include "DifferentialAdhesionGeneralisedLinearSpringForce.hpp" +#include "RandomMotionForce.hpp" + +#include "OnLatticeSimulation.hpp" +#include "CellPopulationAdjacencyMatrixWriter.hpp" + +#include "CaBasedCellPopulation.hpp" +#include "ShovingCaBasedDivisionRule.hpp" +#include "RandomCaSwitchingUpdateRule.hpp" +#include "DifferentialAdhesionCaSwitchingUpdateRule.hpp" + +#include "PottsBasedCellPopulation.hpp" +#include "PottsMeshGenerator.hpp" +#include "VolumeConstraintPottsUpdateRule.hpp" +#include "SurfaceAreaConstraintPottsUpdateRule.hpp" +#include "AdhesionPottsUpdateRule.hpp" +#include "DifferentialAdhesionPottsUpdateRule.hpp" + +#include "CellIdWriter.hpp" +#include "CellMutationStatesWriter.hpp" #include "PetscSetupAndFinalize.hpp" +/* + * This is where you can set parameters to be used in all the simulations. + * + * The first block (commented out) are the original parameter values. + * The second block are parameters for a much shorter simulation, and are used for continuous testing with Chaste. + */ + +//static const double M_TIME_TO_STEADY_STATE = 10; //10 +//static const double M_TIME_FOR_SIMULATION = 100; //100 +//static const double M_NUM_CELLS_ACROSS = 20; //20 // this ^2 cells +//static const double M_CELL_FLUCTUATION = 1.0; +static const double M_TIME_TO_STEADY_STATE = 1.0; //10 +static const double M_TIME_FOR_SIMULATION = 2.0; //100 +static const double M_NUM_CELLS_ACROSS = 20; //20 // this ^2 cells +static const double M_CELL_FLUCTUATION = 1.0; -class Test05CellSorting: public AbstractCellBasedWithTimingsTestSuite +static const std::string M_HEAD_FOLDER = "OpenVT/Test02CellSorting"; + +class Test05CellSorting : public AbstractCellBasedWithTimingsTestSuite { private: + /* + * This is a helper method to randomly label cells add is used in all simulations. + */ + + void RandomlyLabelCells(std::list& rCells, boost::shared_ptr pLabel, double labelledRatio) + { + for (std::list::iterator cell_iter = rCells.begin(); + cell_iter != rCells.end(); + ++cell_iter) + { + if (RandomNumberGenerator::Instance()->ranf() < labelledRatio) + { + (*cell_iter)->AddCellProperty(pLabel); + } + } + } + public: /* + * == CP == + * + * Simulate a population of cells exhibiting cell sorting using the + * Cellular Potts model. */ - void TestNotMuch() - { - TS_ASSERT_THROWS_NOWT(); + void TestPottsMonolayerCellSorting() + { + std::string output_directory = M_HEAD_FOLDER + "/Potts"; + + // Create a simple 2D PottsMesh + unsigned element_size = 4; + unsigned domain_size = M_NUM_CELLS_ACROSS * element_size * 3; // Three times the initial domain size + PottsMeshGenerator<2> generator(domain_size, M_NUM_CELLS_ACROSS, element_size, domain_size, M_NUM_CELLS_ACROSS, element_size); + boost::shared_ptr > p_mesh = generator.GetMesh(); + + // Create cells + std::vector cells; + MAKE_PTR(DifferentiatedCellProliferativeType, p_differentiated_type); + CellsGenerator cells_generator; + cells_generator.GenerateBasicRandom(cells, p_mesh->GetNumElements(), p_differentiated_type); + + // Create cell population + PottsBasedCellPopulation<2> cell_population(*p_mesh, cells); + cell_population.AddCellWriter(); + cell_population.AddCellWriter(); + cell_population.AddPopulationWriter(); + cell_population.AddPopulationWriter(); + + // Set the Temperature + cell_population.SetTemperature(0.2); //Default is 0.1 + + // Set up cell-based simulation and output directory + OnLatticeSimulation<2> simulator(cell_population); + simulator.SetOutputDirectory(output_directory); + + // Set time step and end time for simulation + simulator.SetDt(0.01); + simulator.SetSamplingTimestepMultiple(100); + simulator.SetEndTime(M_TIME_TO_STEADY_STATE); + + // Create update rules and pass to the simulation + MAKE_PTR(VolumeConstraintPottsUpdateRule<2>, p_volume_constraint_update_rule); + p_volume_constraint_update_rule->SetMatureCellTargetVolume(16); // i.e 4x4 cells + p_volume_constraint_update_rule->SetDeformationEnergyParameter(0.1); + simulator.AddUpdateRule(p_volume_constraint_update_rule); + + MAKE_PTR(SurfaceAreaConstraintPottsUpdateRule<2>, p_surface_constraint_update_rule); + p_surface_constraint_update_rule->SetMatureCellTargetSurfaceArea(16); // i.e 4x4 cells + p_surface_constraint_update_rule->SetDeformationEnergyParameter(0.01);//0.01 + simulator.AddUpdateRule(p_surface_constraint_update_rule); + + MAKE_PTR(DifferentialAdhesionPottsUpdateRule<2>, p_differential_adhesion_update_rule); + p_differential_adhesion_update_rule->SetLabelledCellLabelledCellAdhesionEnergyParameter(0.1); + p_differential_adhesion_update_rule->SetLabelledCellCellAdhesionEnergyParameter(0.5); // 1.0 + p_differential_adhesion_update_rule->SetCellCellAdhesionEnergyParameter(0.1); //0.1 + p_differential_adhesion_update_rule->SetCellBoundaryAdhesionEnergyParameter(0.2); // 1.0 + p_differential_adhesion_update_rule->SetLabelledCellBoundaryAdhesionEnergyParameter(1.0); // 2.0 + simulator.AddUpdateRule(p_differential_adhesion_update_rule); + + // Run simulation + simulator.Solve(); + + // Now label some cells + boost::shared_ptr p_state(CellPropertyRegistry::Instance()->Get()); + RandomlyLabelCells(simulator.rGetCellPopulation().rGetCells(), p_state, 0.5); + + // Adjust Parameters + dynamic_cast *>(&(simulator.rGetCellPopulation()))->SetTemperature(0.2*M_CELL_FLUCTUATION); + + // Run simulation + simulator.SetEndTime(M_TIME_TO_STEADY_STATE + M_TIME_FOR_SIMULATION); + simulator.Solve(); + + // Check that the same number of cells + TS_ASSERT_EQUALS(simulator.rGetCellPopulation().GetNumRealCells(), M_NUM_CELLS_ACROSS*M_NUM_CELLS_ACROSS); + + // Test no births or deaths + TS_ASSERT_EQUALS(simulator.GetNumBirths(), 0u); + TS_ASSERT_EQUALS(simulator.GetNumDeaths(), 0u); } + + /* + * == OS == + * + * Simulate a population of cells exhibiting cell sorting using the + * Overlapping Sphere model. + */ + void TestNodeBasedMonolayerCellSorting() + { + std::string output_directory = M_HEAD_FOLDER + "/Node"; + + // Create a simple mesh + HoneycombMeshGenerator generator(M_NUM_CELLS_ACROSS, M_NUM_CELLS_ACROSS, 0); + boost::shared_ptr > p_generating_mesh = generator.GetMesh(); + + //Extended to allow sorting for longer distances + double cut_off_length = 2.5; + + // Convert this to a NodesOnlyMesh + NodesOnlyMesh<2> mesh; + mesh.ConstructNodesWithoutMesh(*p_generating_mesh, cut_off_length); + + // Set up cells, one for each Node + std::vector cells; + MAKE_PTR(DifferentiatedCellProliferativeType, p_differentiated_type); + CellsGenerator cells_generator; + cells_generator.GenerateBasicRandom(cells, mesh.GetNumNodes(), p_differentiated_type); + + // Create cell population + NodeBasedCellPopulation<2> cell_population(mesh, cells); + + // Set population to output all data to results files + cell_population.AddCellWriter(); + cell_population.AddCellWriter(); + cell_population.AddPopulationWriter(); + cell_population.AddPopulationWriter(); + + // Set up cell-based simulation and output directory + OffLatticeSimulation<2> simulator(cell_population); + simulator.SetOutputDirectory(output_directory); + + // Set time step and end time for simulation + simulator.SetDt(1.0/200.0); + simulator.SetSamplingTimestepMultiple(200); + simulator.SetEndTime(M_TIME_TO_STEADY_STATE); + + // Create a force law and pass it to the simulation + MAKE_PTR(DifferentialAdhesionGeneralisedLinearSpringForce<2>, p_differential_adhesion_force); + p_differential_adhesion_force->SetMeinekeSpringStiffness(50.0); + p_differential_adhesion_force->SetHomotypicLabelledSpringConstantMultiplier(1.0); + p_differential_adhesion_force->SetHeterotypicSpringConstantMultiplier(0.1); + p_differential_adhesion_force->SetCutOffLength(cut_off_length); + simulator.AddForce(p_differential_adhesion_force); + + // Add some noise to avoid local minimum + MAKE_PTR(RandomMotionForce<2>, p_random_force); + p_random_force->SetMovementParameter(0.05); //0.1 causes dissasociation, 0.001 is not enough + simulator.AddForce(p_random_force); + + // Run simulation + simulator.Solve(); + + // Now label some cells + boost::shared_ptr p_state(CellPropertyRegistry::Instance()->Get()); + RandomlyLabelCells(simulator.rGetCellPopulation().rGetCells(), p_state, 0.5); + + // Adjust parameters + p_random_force->SetMovementParameter(0.05*M_CELL_FLUCTUATION); //0.1 causes dissasociation + + + // Run simulation + simulator.SetEndTime(M_TIME_TO_STEADY_STATE + M_TIME_FOR_SIMULATION); + simulator.Solve(); + + // Check that the same number of cells + TS_ASSERT_EQUALS(simulator.rGetCellPopulation().GetNumRealCells(), M_NUM_CELLS_ACROSS*M_NUM_CELLS_ACROSS); + + // Test no births or deaths + TS_ASSERT_EQUALS(simulator.GetNumBirths(), 0u); + TS_ASSERT_EQUALS(simulator.GetNumDeaths(), 0u); + } + /* + * == VT == + * + * Simulate a population of cells exhibiting cell sorting using the + * Voronoi tesselation model. + */ + void TestMeshBasedWithGhostsMonolayerCellSorting() + { + std::string output_directory = M_HEAD_FOLDER + "/Mesh"; + + // Create a simple mesh + unsigned num_ghosts = 20; + HoneycombMeshGenerator generator(M_NUM_CELLS_ACROSS, M_NUM_CELLS_ACROSS, num_ghosts); + boost::shared_ptr > p_mesh = generator.GetMesh(); + + // Set up cells, one for each non ghost Node + std::vector location_indices = generator.GetCellLocationIndices(); + std::vector cells; + MAKE_PTR(DifferentiatedCellProliferativeType, p_differentiated_type); + CellsGenerator cells_generator; + cells_generator.GenerateBasicRandom(cells, location_indices.size(), p_differentiated_type); + + // Create cell population + MeshBasedCellPopulationWithGhostNodes<2> cell_population(*p_mesh, cells, location_indices); + cell_population.AddPopulationWriter(); + + // Set population to output all data to results files + cell_population.AddCellWriter(); + cell_population.AddCellWriter(); + cell_population.AddPopulationWriter(); + + // Set up cell-based simulation and output directory + OffLatticeSimulation<2> simulator(cell_population); + simulator.SetOutputDirectory(output_directory); + + // Set time step and end time for simulation + simulator.SetDt(1.0/200.0); + simulator.SetSamplingTimestepMultiple(200); + simulator.SetEndTime(M_TIME_TO_STEADY_STATE); + + // Create a force law and pass it to the simulation + MAKE_PTR(DifferentialAdhesionGeneralisedLinearSpringForce<2>, p_differential_adhesion_force); + p_differential_adhesion_force->SetMeinekeSpringStiffness(50.0); + p_differential_adhesion_force->SetHomotypicLabelledSpringConstantMultiplier(1.0); + p_differential_adhesion_force->SetHeterotypicSpringConstantMultiplier(0.1); + simulator.AddForce(p_differential_adhesion_force); + + // Add some noise to avoid local minimum + MAKE_PTR(RandomMotionForce<2>, p_random_force); + p_random_force->SetMovementParameter(0.1); + simulator.AddForce(p_random_force); + + // Run simulation + simulator.Solve(); + + // Now label some cells + boost::shared_ptr p_state(CellPropertyRegistry::Instance()->Get()); + RandomlyLabelCells(simulator.rGetCellPopulation().rGetCells(), p_state, 0.5); + + // Adjust parameters + p_random_force->SetMovementParameter(0.1*M_CELL_FLUCTUATION); //0.1 causes dissasociation + + // Run simulation + simulator.SetEndTime(M_TIME_TO_STEADY_STATE + M_TIME_FOR_SIMULATION); + simulator.Solve(); + + // Check that the same number of cells + TS_ASSERT_EQUALS(simulator.rGetCellPopulation().GetNumRealCells(), M_NUM_CELLS_ACROSS*M_NUM_CELLS_ACROSS); + + // Test no births or deaths + TS_ASSERT_EQUALS(simulator.GetNumBirths(), 0u); + TS_ASSERT_EQUALS(simulator.GetNumDeaths(), 0u); + } + + /* + * == VM == + * + * Simulate a population of cells exhibiting cell sorting using the + * Cell Vertex model. + */ + void TestVertexMonolayerCellSorting() + { + std::string output_directory = M_HEAD_FOLDER + "/Vertex"; + + // Create a simple 2D MutableVertexMesh + HoneycombVertexMeshGenerator generator(M_NUM_CELLS_ACROSS, M_NUM_CELLS_ACROSS); + boost::shared_ptr > p_mesh = generator.GetMesh(); + p_mesh->SetCellRearrangementThreshold(0.1); + + // Slows things down but can use a larger timestep and diffusion forces + //p_mesh->SetCheckForInternalIntersections(true); + + // Set up cells, one for each VertexElement + std::vector cells; + boost::shared_ptr p_cell_type(CellPropertyRegistry::Instance()->Get()); + CellsGenerator cells_generator; + cells_generator.GenerateBasicRandom(cells, p_mesh->GetNumElements(), p_cell_type); + + for (unsigned i=0; iGetCellData()->SetItem("target area", 1.0); + } + + // Create cell population + VertexBasedCellPopulation<2> cell_population(*p_mesh, cells); + + // Set population to output all data to results files + cell_population.AddCellWriter(); + cell_population.AddCellWriter(); + cell_population.AddPopulationWriter(); + cell_population.AddPopulationWriter(); + + // Set up cell-based simulation and output directory + OffLatticeSimulation<2> simulator(cell_population); + simulator.SetOutputDirectory(output_directory); + + // Set time step and end time for simulation + simulator.SetDt(1.0/200.0); + simulator.SetSamplingTimestepMultiple(200); + simulator.SetEndTime(M_TIME_TO_STEADY_STATE); + + // Set up force law and pass it to the simulation + MAKE_PTR(NagaiHondaDifferentialAdhesionForce<2>, p_force); + p_force->SetNagaiHondaDeformationEnergyParameter(50.0); + p_force->SetNagaiHondaMembraneSurfaceEnergyParameter(1.0); + p_force->SetNagaiHondaCellCellAdhesionEnergyParameter(1.0); + p_force->SetNagaiHondaLabelledCellCellAdhesionEnergyParameter(2.0); + p_force->SetNagaiHondaLabelledCellLabelledCellAdhesionEnergyParameter(1.0); + p_force->SetNagaiHondaCellBoundaryAdhesionEnergyParameter(10.0); + p_force->SetNagaiHondaLabelledCellBoundaryAdhesionEnergyParameter(20.0); + simulator.AddForce(p_force); + + // Add some noise to avoid local minimum + MAKE_PTR(RandomMotionForce<2>, p_random_force); + p_random_force->SetMovementParameter(0.1); + simulator.AddForce(p_random_force); + + // Run simulation + simulator.Solve(); + + // Now label some cells + boost::shared_ptr p_state(CellPropertyRegistry::Instance()->Get()); + RandomlyLabelCells(simulator.rGetCellPopulation().rGetCells(), p_state, 0.5); + + // Adjust parameters + p_random_force->SetMovementParameter(0.1*M_CELL_FLUCTUATION); + + // Run simulation + simulator.SetEndTime(M_TIME_TO_STEADY_STATE + M_TIME_FOR_SIMULATION); + simulator.Solve(); + + + // Check that the same number of cells + TS_ASSERT_EQUALS(simulator.rGetCellPopulation().GetNumRealCells(), M_NUM_CELLS_ACROSS*M_NUM_CELLS_ACROSS); + + // Test no births or deaths + TS_ASSERT_EQUALS(simulator.GetNumBirths(), 0u); + TS_ASSERT_EQUALS(simulator.GetNumDeaths(), 0u); + } + }; #endif /* TEST05CELLSORTING_HPP_ */