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

Simulation of fission in epithelial layer

Introduction

In this test we show how Chaste can be used to simulate a buckling layer of epithelial cells. Details of the computational model can be found in Langlands et al (2016) “Paneth cell-rich regions separated by a cluster of Lgr5+ cells initiate fission in the intestinal stem cell niche”.

Including header files

We begin by including the necessary header files. The first ones are common to all cell_based Chaste simulations

#include <cxxtest/TestSuite.h> //Needed for all test files
#include "CellBasedSimulationArchiver.hpp" //Needed if we would like to save/load simulations
#include "AbstractCellBasedTestSuite.hpp" //Needed for cell-based tests: times simulations, generates random numbers and has cell properties
#include "CheckpointArchiveTypes.hpp" //Needed if we use GetIdentifier() method (which we do)
#include "SmartPointers.hpp" //Enables macros to save typing

The next set of classes are needed specifically for the simulation, which can be found in the core code.

#include "HoneycombMeshGenerator.hpp" //Generates mesh
#include "OffLatticeSimulation.hpp" //Simulates the evolution of the population
#include "MeshBasedCellPopulationWithGhostNodes.hpp"
#include "VoronoiDataWriter.hpp" //Allows us to visualise output in Paraview
#include "DifferentiatedCellProliferativeType.hpp" //Stops cells from proliferating
#include "TransitCellProliferativeType.hpp"
#include "FakePetscSetup.hpp" //Forbids tests running in parallel

This set of classes are not part of the core code and were made for this model.

#include "StochasticTargetProportionBasedCellCycleModel.hpp" //Asymmetric-division-based cell cycle model
#include "PanethCellMutationState.hpp" //Mutation class that defines Paneth cells

#include "EpithelialLayerLinearSpringForce.hpp" //Spring force law to account for different cell type pairs
#include "EpithelialLayerBasementMembraneForce.hpp" //Basement membrane force, as defined in Dunn et al. (2012)
#include "EpithelialLayerAnoikisCellKiller.hpp" //Cell killer to remove proliferative cells that have fallen into the lumen

#include "EpithelialLayerDataTrackingModifier.hpp" //Modifier for all the necessary data

#include "FixedRegionPlaneBoundaryCondition.hpp" //Boundary condition that fixes cells past a given plane

Define the Chaste simulation as a test class. This is how all simulations in Chaste are defined.

class TestCryptFissionLiteratePaper : public AbstractCellBasedTestSuite
{
public:
	void TestEpithelialLayerUndergoingFission() throw(Exception)
	{

We first set all the simulation parameters.

		//Simulation time parameters
		double dt = 0.005; //Set dt
		double end_time = 100.0; //Set end time
		double sampling_timestep = 0.5/dt; //Set sampling timestep

		//Set all the spring stiffness variables
		double epithelial_epithelial_stiffness = 15.0; //Epithelial-epithelial spring connections
		double epithelial_nonepithelial_stiffness = 15.0; //Epithelial-non-epithelial spring connections
		double nonepithelial_nonepithelial_stiffness = 15.0; //Non-epithelial-non-epithelial spring connections

		//Set the stiffness ratio for Paneth cells to stem cells. This is the
		double stiffness_ratio = 4.5;

Set the target proportion for stem cells. Note that the target proportion and the actual proportion of stem cells will be different, due to cell death and mechanical heterogeneities.

		double target_proportion = 0.3; //This corresponds to the 20% stem cell (80% Paneth cell) case
		// double target_proportion = 0.8; //This corresponds to the 60% stem cell (40% Paneth cell) case.

		//Set the BM force parameters
		double bm_force = 10.0; //Set the basement membrane stiffness
		double target_curvature = 0.2; //Set the target curvature, i.e. how circular the layer wants to be

Set the domain of the model.

		unsigned cells_across = 24; //Desired width + a few more layers, to avoid the box collapsing
		unsigned cells_up = 27; //Since height of each cell is 0.5*sqrt(3), we need to specify the no. of cells such that #cells*0.5*sqrt(3) = desired height
		unsigned ghosts = 4; //Define a sufficient layer of ghost nodes to avoid infinite tessellations and hence excessively large forces

		//Translate mesh 1.5 units left and 1.5 units down, so that we have a sufficient layer of fixed cells around the boundary.
		c_vector<double, 2> translate_left = zero_vector<double>(2);
		translate_left(0) = -1.5;

		c_vector<double, 2> translate_down = zero_vector<double>(2);
		translate_down(1) = -1.5;

Define the initially circular lumen by centre and radius

		c_vector<double,2> circle_centre;
		circle_centre(0) = 10.5;
		circle_centre(1) = 10.0;

		double circle_radius = 2.5; //Size of hole
		assert(circle_radius > 0); //Just in case someone does something crazy.

		double ring_radius = circle_radius + 2.0; //Radius of the ring of cells. This isn't the actual radius, just has to be large enough for later
		assert((ring_radius <= cells_across)&&(ring_radius <= cells_up)); //Again, just in case.

Generate the initial mesh of cells.

		HoneycombMeshGenerator generator(cells_across, cells_up, ghosts);
		MutableMesh<2,2>* p_mesh = generator.GetMesh();

		//Translate mesh appropriately
		p_mesh->Translate(translate_left);
		p_mesh->Translate(translate_down);

Define the lumen as an inner region of ghost nodes.

		std::vector<unsigned> initial_real_indices = generator.GetCellLocationIndices(); //Obtain the locations of real nodes

		std::vector<unsigned> real_indices; //Vector used to define the locations of non-ghost nodes

		//Sweep over the initial real indices
		for (unsigned i = 0; i < initial_real_indices.size(); i++)
		{
			unsigned cell_index = initial_real_indices[i];
			double x = p_mesh->GetNode(cell_index)->rGetLocation()[0];
			double y = p_mesh->GetNode(cell_index)->rGetLocation()[1];

			// If the location of the node falls inside the defined lumen region, then it becomes a ghost node.
			if (pow(x-circle_centre[0],2) + pow(y-circle_centre[1],2) > pow(circle_radius,2))
			{
				real_indices.push_back(cell_index);
			}
		}

Define cell types: non-epithelial cells are differentiated, all proliferative epithelial cells will be transit cells (to define cell cycle duration). In addition, stem cells are assigned a wildtype mutation state, while Paneth cells are assigned a Paneth cell mutation state.

		boost::shared_ptr<AbstractCellProperty> p_diff_type = CellPropertyRegistry::Instance()->Get<DifferentiatedCellProliferativeType>();
		boost::shared_ptr<AbstractCellProperty> p_stem_type = CellPropertyRegistry::Instance()->Get<TransitCellProliferativeType>();
		boost::shared_ptr<AbstractCellProperty> p_paneth_state = CellPropertyRegistry::Instance()->Get<PanethCellMutationState>();
		boost::shared_ptr<AbstractCellProperty> p_state = CellPropertyRegistry::Instance()->Get<WildTypeCellMutationState>();

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

Create asymmetric-division-based cell cycle for each cell. However, we initially set all cells to be non-epithelial cells before defining our layer of epithelial cells.

		for (unsigned i = 0; i<real_indices.size(); i++)
		{
			//Set cell cycle
			StochasticTargetProportionBasedCellCycleModel* p_cycle_model = new StochasticTargetProportionBasedCellCycleModel();
			p_cycle_model->SetTargetProportion(target_proportion); //Set the division parameter
			p_cycle_model->SetDimension(2);

			//To avoid a 'pulsing' behaviour with birth events, we set each cell's initial age to be
			// ~U(-12, 0) in the past, as each cell cycle duration is U(11, 13).
			double birth_time = 12.0*RandomNumberGenerator::Instance()->ranf();
			p_cycle_model->SetBirthTime(-birth_time);

			CellPtr p_cell(new Cell(p_state, p_cycle_model));
			p_cell->SetCellProliferativeType(p_diff_type); //Set the cell to be differentiated and hence non-epithelial
			p_cell->InitialiseCellCycleModel();

			cells.push_back(p_cell);
		}

		//Create cell population
		MeshBasedCellPopulationWithGhostNodes<2> cell_population(*p_mesh, cells, real_indices);

		//Create layer of proliferative cells
		for (unsigned i = 0; i < real_indices.size(); i++)
		{
			unsigned cell_index = real_indices[i];
			CellPtr cell_iter = cell_population.GetCellUsingLocationIndex(cell_index);
			double x = cell_population.GetLocationOfCellCentre(cell_iter)[0];
			double y = cell_population.GetLocationOfCellCentre(cell_iter)[1];

We ‘un-differentiate’ any cells adjacent to the lumen into epithelial cells. We only consider cells within the pre-defined ring radius. Not only does this narrow down our search, but it also means we won’t accidentally consider any cells on the outside and turn them into epithelial cells.

			if (pow(x-circle_centre[0],2) + pow(y-circle_centre[1],2) <= pow(ring_radius,2))
			{
				Node<2>* p_node = cell_population.GetNodeCorrespondingToCell(cell_iter);
				unsigned node_index = p_node->GetIndex();

				//Iterate over all possible neighbours of the node
				for (Node<2>::ContainingElementIterator iter = p_node->ContainingElementsBegin();
						iter != p_node->ContainingElementsEnd();
						++iter)
				{
					bool element_contains_ghost_nodes = false;

					// Get a pointer to the element
					Element<2,2>* p_element = cell_population.rGetMesh().GetElement(*iter);

					// Check whether it's triangulation contains a ghost node
					for (unsigned local_index=0; local_index<3; local_index++)
					{
						unsigned nodeGlobalIndex = p_element->GetNodeGlobalIndex(local_index);

						if (cell_population.IsGhostNode(nodeGlobalIndex) == true)
						{
							element_contains_ghost_nodes = true;
							break; 				// This should break out of the inner for loop
						}
					}

					//If a cell has a ghost node as a neighbour, we make it an epithelial cells.
					if(element_contains_ghost_nodes)
					{
						cell_iter->SetCellProliferativeType(p_stem_type);
					}
				}
			}
		}

Iterate again and check that proliferative cells are also attached to non-epithelial cells. If they are not, remove them from the simulation.

		for (unsigned i = 0; i < real_indices.size(); i++)
		{
			unsigned cell_index = real_indices[i];
			CellPtr cell_iter = cell_population.GetCellUsingLocationIndex(cell_index);
			double x = cell_population.GetLocationOfCellCentre(cell_iter)[0];
			double y = cell_population.GetLocationOfCellCentre(cell_iter)[1];

			//Only consider this inside the pre-defined ring to narrow down our search

			if (pow(x-circle_centre[0],2) + pow(y-circle_centre[1],2) <= pow(ring_radius,2))
			{
				Node<2>* p_node = cell_population.GetNodeCorrespondingToCell(cell_iter);
				unsigned node_index = p_node->GetIndex();

				//Only iterate over the initial layer of transit cells
				if (cell_iter->GetCellProliferativeType()->IsType<DifferentiatedCellProliferativeType>() == false)
				{
					bool element_contains_gel_nodes = false;

					//Iterate over elements (triangles) containing the node
					for (Node<2>::ContainingElementIterator iter = p_node->ContainingElementsBegin();
							iter != p_node->ContainingElementsEnd();
							++iter)
					{
						// Get a pointer to the element (triangle)
						Element<2,2>* p_element = cell_population.rGetMesh().GetElement(*iter);

						// Check if its triangulation contains a gel node
						for (unsigned local_index=0; local_index<3; local_index++)
						{
							unsigned nodeGlobalIndex = p_element->GetNodeGlobalIndex(local_index);
							bool is_ghost_node = cell_population.IsGhostNode(nodeGlobalIndex);

							if (is_ghost_node == false) //Make sure we're not dealing with ghost nodes (otherwise this stuff will fail)
							{
								CellPtr p_local_cell = cell_population.GetCellUsingLocationIndex(nodeGlobalIndex);
								if (p_local_cell->GetCellProliferativeType()->IsType<DifferentiatedCellProliferativeType>()==true)
								{
									element_contains_gel_nodes = true;
									break; 				// This should break out of the inner for loop
								}
							}
						}
					}

					if(element_contains_gel_nodes == false)
					{
						cell_iter->Kill();
					}
				}
			}
		}

Randomly assign cells in the layer to be Paneth cells.

		std::vector<unsigned> cells_in_layer; //Initialise vector

		//Obtain the proliferative cells
		for (AbstractCellPopulation<2>::Iterator cell_iter = cell_population.Begin();
				cell_iter != cell_population.End();
				++cell_iter)
		{
			unsigned node_index = cell_population.GetLocationIndexUsingCell(*cell_iter);

			//If the cell is an epithelial cell
			if (!cell_iter->GetCellProliferativeType()->IsType<DifferentiatedCellProliferativeType>() )
			{
				cells_in_layer.push_back(node_index); //Add the angle and node index
			}
		}

For each cell in the ring, we draw a random number and assign cells to be stem cells with a probability equal to the target proportion, as defined above.

		for (unsigned i = 0; i < cells_in_layer.size(); i++)
		{
			unsigned node_index = cells_in_layer[i];

			CellPtr cell = cell_population.GetCellUsingLocationIndex(node_index);

			//Randomly generate number
			double random_number = RandomNumberGenerator::Instance()->ranf();

			if(random_number >= target_proportion) //Assign cells to be Paneth with 1 - target_proportion
			{
				cell->SetMutationState(p_paneth_state);
			}
		}

		//Allow output in Paraview, a program that can be used to visualise Chaste simulations
		cell_population.AddPopulationWriter<VoronoiDataWriter>();

Define the simulation class.

		OffLatticeSimulation<2> simulator(cell_population);

		//Set output directory
		simulator.SetOutputDirectory("CryptFissionLiteratePaper");

		simulator.SetDt(dt); //Set the timestep dt for force volution
		simulator.SetSamplingTimestepMultiple(sampling_timestep); //Set the sampling timestep multiple for animations
		simulator.SetEndTime(end_time); //Set the number of hours to run the simulation to

We add a modifier class to track relevant cell population numbers and shape measurements.

		MAKE_PTR(EpithelialLayerDataTrackingModifier<2>, p_data_tracking_modifier);
		simulator.AddSimulationModifier(p_data_tracking_modifier);

Add linear spring force which has different spring stiffness constants, depending on the pair of cells it is connecting.

		MAKE_PTR(EpithelialLayerLinearSpringForce<2>, p_spring_force);
		p_spring_force->SetCutOffLength(1.5);
		//Set the spring stiffnesses
		p_spring_force->SetEpithelialEpithelialSpringStiffness(epithelial_epithelial_stiffness);
		p_spring_force->SetEpithelialNonepithelialSpringStiffness(epithelial_nonepithelial_stiffness);
		p_spring_force->SetNonepithelialNonepithelialSpringStiffness(nonepithelial_nonepithelial_stiffness);
		p_spring_force->SetPanethCellStiffnessRatio(stiffness_ratio);
		simulator.AddForce(p_spring_force);

Add the basement membrane force.

		MAKE_PTR(EpithelialLayerBasementMembraneForce, p_bm_force);
		p_bm_force->SetBasementMembraneParameter(bm_force); //Equivalent to beta in SJD's papers
		p_bm_force->SetTargetCurvature(target_curvature); //This is equivalent to 1/R in SJD's papers
		simulator.AddForce(p_bm_force);

Add an anoikis-based cell killer.

		MAKE_PTR_ARGS(EpithelialLayerAnoikisCellKiller, p_anoikis_killer, (&cell_population));
		simulator.AddCellKiller(p_anoikis_killer);

We fix all cells outside of the 20 x 20 box.

		c_vector<double,2> point = zero_vector<double>(2);
		c_vector<double,2> normal = zero_vector<double>(2);

		//Fix cells in the region x < 0
		normal(0) = -1.0;
		MAKE_PTR_ARGS(FixedRegionPlaneBoundaryCondition<2>, p_bc1, (&cell_population, point, normal));
		simulator.AddCellPopulationBoundaryCondition(p_bc1);

		//Fix cells in the region x > 20
		point(0) = 20.0;
		normal(0) = 1.0;
		MAKE_PTR_ARGS(FixedRegionPlaneBoundaryCondition<2>, p_bc2, (&cell_population, point, normal));
		simulator.AddCellPopulationBoundaryCondition(p_bc2);

		//Fix cells in the region y < 0
		point(0) = 0.0;
		point(1) = 0.0;
		normal(0) = 0.0;
		normal(1) = -1.0;
		MAKE_PTR_ARGS(FixedRegionPlaneBoundaryCondition<2>, p_bc3, (&cell_population, point, normal));
		simulator.AddCellPopulationBoundaryCondition(p_bc3);

		//Fix cells in the region y > 20
		point(1) = 20.0;
		normal(1) = 1.0;
		MAKE_PTR_ARGS(FixedRegionPlaneBoundaryCondition<2>, p_bc4, (&cell_population, point, normal));
		simulator.AddCellPopulationBoundaryCondition(p_bc4);

Run the simulation.

		simulator.Solve();
	}

To visualize the results, open a new terminal and

cd

to the Chaste directory, then

cd

to

anim

. Then do:

java Visualize2dCentreCells /tmp/$USER/testoutput/CryptFissionLiteratePaper/results_from_time_0

. You may have to do:

javac Visualize2dCentreCells.java

beforehand to create the java executable. You should also select the axes equal option.

};

Code

The full code is given below

File name TestCryptFissionLiteratePaper.hpp

#include <cxxtest/TestSuite.h> //Needed for all test files
#include "CellBasedSimulationArchiver.hpp" //Needed if we would like to save/load simulations
#include "AbstractCellBasedTestSuite.hpp" //Needed for cell-based tests: times simulations, generates random numbers and has cell properties
#include "CheckpointArchiveTypes.hpp" //Needed if we use GetIdentifier() method (which we do)
#include "SmartPointers.hpp" //Enables macros to save typing

#include "HoneycombMeshGenerator.hpp" //Generates mesh
#include "OffLatticeSimulation.hpp" //Simulates the evolution of the population
#include "MeshBasedCellPopulationWithGhostNodes.hpp"
#include "VoronoiDataWriter.hpp" //Allows us to visualise output in Paraview
#include "DifferentiatedCellProliferativeType.hpp" //Stops cells from proliferating
#include "TransitCellProliferativeType.hpp"
#include "FakePetscSetup.hpp" //Forbids tests running in parallel

#include "StochasticTargetProportionBasedCellCycleModel.hpp" //Asymmetric-division-based cell cycle model
#include "PanethCellMutationState.hpp" //Mutation class that defines Paneth cells

#include "EpithelialLayerLinearSpringForce.hpp" //Spring force law to account for different cell type pairs
#include "EpithelialLayerBasementMembraneForce.hpp" //Basement membrane force, as defined in Dunn et al. (2012)
#include "EpithelialLayerAnoikisCellKiller.hpp" //Cell killer to remove proliferative cells that have fallen into the lumen

#include "EpithelialLayerDataTrackingModifier.hpp" //Modifier for all the necessary data

#include "FixedRegionPlaneBoundaryCondition.hpp" //Boundary condition that fixes cells past a given plane

class TestCryptFissionLiteratePaper : public AbstractCellBasedTestSuite
{
public:
	void TestEpithelialLayerUndergoingFission() throw(Exception)
	{
		//Simulation time parameters
		double dt = 0.005; //Set dt
		double end_time = 100.0; //Set end time
		double sampling_timestep = 0.5/dt; //Set sampling timestep

		//Set all the spring stiffness variables
		double epithelial_epithelial_stiffness = 15.0; //Epithelial-epithelial spring connections
		double epithelial_nonepithelial_stiffness = 15.0; //Epithelial-non-epithelial spring connections
		double nonepithelial_nonepithelial_stiffness = 15.0; //Non-epithelial-non-epithelial spring connections

		//Set the stiffness ratio for Paneth cells to stem cells. This is the
		double stiffness_ratio = 4.5;

		double target_proportion = 0.3; //This corresponds to the 20% stem cell (80% Paneth cell) case
		// double target_proportion = 0.8; //This corresponds to the 60% stem cell (40% Paneth cell) case.

		//Set the BM force parameters
		double bm_force = 10.0; //Set the basement membrane stiffness
		double target_curvature = 0.2; //Set the target curvature, i.e. how circular the layer wants to be

		unsigned cells_across = 24; //Desired width + a few more layers, to avoid the box collapsing
		unsigned cells_up = 27; //Since height of each cell is 0.5*sqrt(3), we need to specify the no. of cells such that #cells*0.5*sqrt(3) = desired height
		unsigned ghosts = 4; //Define a sufficient layer of ghost nodes to avoid infinite tessellations and hence excessively large forces

		//Translate mesh 1.5 units left and 1.5 units down, so that we have a sufficient layer of fixed cells around the boundary.
		c_vector<double, 2> translate_left = zero_vector<double>(2);
		translate_left(0) = -1.5;

		c_vector<double, 2> translate_down = zero_vector<double>(2);
		translate_down(1) = -1.5;

		c_vector<double,2> circle_centre;
		circle_centre(0) = 10.5;
		circle_centre(1) = 10.0;

		double circle_radius = 2.5; //Size of hole
		assert(circle_radius > 0); //Just in case someone does something crazy.

		double ring_radius = circle_radius + 2.0; //Radius of the ring of cells. This isn't the actual radius, just has to be large enough for later
		assert((ring_radius <= cells_across)&&(ring_radius <= cells_up)); //Again, just in case.

		HoneycombMeshGenerator generator(cells_across, cells_up, ghosts);
		MutableMesh<2,2>* p_mesh = generator.GetMesh();

		//Translate mesh appropriately
		p_mesh->Translate(translate_left);
		p_mesh->Translate(translate_down);

		std::vector<unsigned> initial_real_indices = generator.GetCellLocationIndices(); //Obtain the locations of real nodes

		std::vector<unsigned> real_indices; //Vector used to define the locations of non-ghost nodes

		//Sweep over the initial real indices
		for (unsigned i = 0; i < initial_real_indices.size(); i++)
		{
			unsigned cell_index = initial_real_indices[i];
			double x = p_mesh->GetNode(cell_index)->rGetLocation()[0];
			double y = p_mesh->GetNode(cell_index)->rGetLocation()[1];

			// If the location of the node falls inside the defined lumen region, then it becomes a ghost node.
			if (pow(x-circle_centre[0],2) + pow(y-circle_centre[1],2) > pow(circle_radius,2))
			{
				real_indices.push_back(cell_index);
			}
		}

		boost::shared_ptr<AbstractCellProperty> p_diff_type = CellPropertyRegistry::Instance()->Get<DifferentiatedCellProliferativeType>();
		boost::shared_ptr<AbstractCellProperty> p_stem_type = CellPropertyRegistry::Instance()->Get<TransitCellProliferativeType>();
		boost::shared_ptr<AbstractCellProperty> p_paneth_state = CellPropertyRegistry::Instance()->Get<PanethCellMutationState>();
		boost::shared_ptr<AbstractCellProperty> p_state = CellPropertyRegistry::Instance()->Get<WildTypeCellMutationState>();

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

		for (unsigned i = 0; i<real_indices.size(); i++)
		{
			//Set cell cycle
			StochasticTargetProportionBasedCellCycleModel* p_cycle_model = new StochasticTargetProportionBasedCellCycleModel();
			p_cycle_model->SetTargetProportion(target_proportion); //Set the division parameter
			p_cycle_model->SetDimension(2);

			//To avoid a 'pulsing' behaviour with birth events, we set each cell's initial age to be
			// ~U(-12, 0) in the past, as each cell cycle duration is U(11, 13).
			double birth_time = 12.0*RandomNumberGenerator::Instance()->ranf();
			p_cycle_model->SetBirthTime(-birth_time);

			CellPtr p_cell(new Cell(p_state, p_cycle_model));
			p_cell->SetCellProliferativeType(p_diff_type); //Set the cell to be differentiated and hence non-epithelial
			p_cell->InitialiseCellCycleModel();

			cells.push_back(p_cell);
		}

		//Create cell population
		MeshBasedCellPopulationWithGhostNodes<2> cell_population(*p_mesh, cells, real_indices);

		//Create layer of proliferative cells
		for (unsigned i = 0; i < real_indices.size(); i++)
		{
			unsigned cell_index = real_indices[i];
			CellPtr cell_iter = cell_population.GetCellUsingLocationIndex(cell_index);
			double x = cell_population.GetLocationOfCellCentre(cell_iter)[0];
			double y = cell_population.GetLocationOfCellCentre(cell_iter)[1];

			if (pow(x-circle_centre[0],2) + pow(y-circle_centre[1],2) <= pow(ring_radius,2))
			{
				Node<2>* p_node = cell_population.GetNodeCorrespondingToCell(cell_iter);
				unsigned node_index = p_node->GetIndex();

				//Iterate over all possible neighbours of the node
				for (Node<2>::ContainingElementIterator iter = p_node->ContainingElementsBegin();
						iter != p_node->ContainingElementsEnd();
						++iter)
				{
					bool element_contains_ghost_nodes = false;

					// Get a pointer to the element
					Element<2,2>* p_element = cell_population.rGetMesh().GetElement(*iter);

					// Check whether it's triangulation contains a ghost node
					for (unsigned local_index=0; local_index<3; local_index++)
					{
						unsigned nodeGlobalIndex = p_element->GetNodeGlobalIndex(local_index);

						if (cell_population.IsGhostNode(nodeGlobalIndex) == true)
						{
							element_contains_ghost_nodes = true;
							break; 				// This should break out of the inner for loop
						}
					}

					//If a cell has a ghost node as a neighbour, we make it an epithelial cells.
					if(element_contains_ghost_nodes)
					{
						cell_iter->SetCellProliferativeType(p_stem_type);
					}
				}
			}
		}

		for (unsigned i = 0; i < real_indices.size(); i++)
		{
			unsigned cell_index = real_indices[i];
			CellPtr cell_iter = cell_population.GetCellUsingLocationIndex(cell_index);
			double x = cell_population.GetLocationOfCellCentre(cell_iter)[0];
			double y = cell_population.GetLocationOfCellCentre(cell_iter)[1];

			//Only consider this inside the pre-defined ring to narrow down our search

			if (pow(x-circle_centre[0],2) + pow(y-circle_centre[1],2) <= pow(ring_radius,2))
			{
				Node<2>* p_node = cell_population.GetNodeCorrespondingToCell(cell_iter);
				unsigned node_index = p_node->GetIndex();

				//Only iterate over the initial layer of transit cells
				if (cell_iter->GetCellProliferativeType()->IsType<DifferentiatedCellProliferativeType>() == false)
				{
					bool element_contains_gel_nodes = false;

					//Iterate over elements (triangles) containing the node
					for (Node<2>::ContainingElementIterator iter = p_node->ContainingElementsBegin();
							iter != p_node->ContainingElementsEnd();
							++iter)
					{
						// Get a pointer to the element (triangle)
						Element<2,2>* p_element = cell_population.rGetMesh().GetElement(*iter);

						// Check if its triangulation contains a gel node
						for (unsigned local_index=0; local_index<3; local_index++)
						{
							unsigned nodeGlobalIndex = p_element->GetNodeGlobalIndex(local_index);
							bool is_ghost_node = cell_population.IsGhostNode(nodeGlobalIndex);

							if (is_ghost_node == false) //Make sure we're not dealing with ghost nodes (otherwise this stuff will fail)
							{
								CellPtr p_local_cell = cell_population.GetCellUsingLocationIndex(nodeGlobalIndex);
								if (p_local_cell->GetCellProliferativeType()->IsType<DifferentiatedCellProliferativeType>()==true)
								{
									element_contains_gel_nodes = true;
									break; 				// This should break out of the inner for loop
								}
							}
						}
					}

					if(element_contains_gel_nodes == false)
					{
						cell_iter->Kill();
					}
				}
			}
		}

		std::vector<unsigned> cells_in_layer; //Initialise vector

		//Obtain the proliferative cells
		for (AbstractCellPopulation<2>::Iterator cell_iter = cell_population.Begin();
				cell_iter != cell_population.End();
				++cell_iter)
		{
			unsigned node_index = cell_population.GetLocationIndexUsingCell(*cell_iter);

			//If the cell is an epithelial cell
			if (!cell_iter->GetCellProliferativeType()->IsType<DifferentiatedCellProliferativeType>() )
			{
				cells_in_layer.push_back(node_index); //Add the angle and node index
			}
		}

		for (unsigned i = 0; i < cells_in_layer.size(); i++)
		{
			unsigned node_index = cells_in_layer[i];

			CellPtr cell = cell_population.GetCellUsingLocationIndex(node_index);

			//Randomly generate number
			double random_number = RandomNumberGenerator::Instance()->ranf();

			if(random_number >= target_proportion) //Assign cells to be Paneth with 1 - target_proportion
			{
				cell->SetMutationState(p_paneth_state);
			}
		}

		//Allow output in Paraview, a program that can be used to visualise Chaste simulations
		cell_population.AddPopulationWriter<VoronoiDataWriter>();

		OffLatticeSimulation<2> simulator(cell_population);

		//Set output directory
		simulator.SetOutputDirectory("CryptFissionLiteratePaper");

		simulator.SetDt(dt); //Set the timestep dt for force volution
		simulator.SetSamplingTimestepMultiple(sampling_timestep); //Set the sampling timestep multiple for animations
		simulator.SetEndTime(end_time); //Set the number of hours to run the simulation to

		MAKE_PTR(EpithelialLayerDataTrackingModifier<2>, p_data_tracking_modifier);
		simulator.AddSimulationModifier(p_data_tracking_modifier);

		MAKE_PTR(EpithelialLayerLinearSpringForce<2>, p_spring_force);
		p_spring_force->SetCutOffLength(1.5);
		//Set the spring stiffnesses
		p_spring_force->SetEpithelialEpithelialSpringStiffness(epithelial_epithelial_stiffness);
		p_spring_force->SetEpithelialNonepithelialSpringStiffness(epithelial_nonepithelial_stiffness);
		p_spring_force->SetNonepithelialNonepithelialSpringStiffness(nonepithelial_nonepithelial_stiffness);
		p_spring_force->SetPanethCellStiffnessRatio(stiffness_ratio);
		simulator.AddForce(p_spring_force);

		MAKE_PTR(EpithelialLayerBasementMembraneForce, p_bm_force);
		p_bm_force->SetBasementMembraneParameter(bm_force); //Equivalent to beta in SJD's papers
		p_bm_force->SetTargetCurvature(target_curvature); //This is equivalent to 1/R in SJD's papers
		simulator.AddForce(p_bm_force);

		MAKE_PTR_ARGS(EpithelialLayerAnoikisCellKiller, p_anoikis_killer, (&cell_population));
		simulator.AddCellKiller(p_anoikis_killer);

		c_vector<double,2> point = zero_vector<double>(2);
		c_vector<double,2> normal = zero_vector<double>(2);

		//Fix cells in the region x < 0
		normal(0) = -1.0;
		MAKE_PTR_ARGS(FixedRegionPlaneBoundaryCondition<2>, p_bc1, (&cell_population, point, normal));
		simulator.AddCellPopulationBoundaryCondition(p_bc1);

		//Fix cells in the region x > 20
		point(0) = 20.0;
		normal(0) = 1.0;
		MAKE_PTR_ARGS(FixedRegionPlaneBoundaryCondition<2>, p_bc2, (&cell_population, point, normal));
		simulator.AddCellPopulationBoundaryCondition(p_bc2);

		//Fix cells in the region y < 0
		point(0) = 0.0;
		point(1) = 0.0;
		normal(0) = 0.0;
		normal(1) = -1.0;
		MAKE_PTR_ARGS(FixedRegionPlaneBoundaryCondition<2>, p_bc3, (&cell_population, point, normal));
		simulator.AddCellPopulationBoundaryCondition(p_bc3);

		//Fix cells in the region y > 20
		point(1) = 20.0;
		normal(1) = 1.0;
		MAKE_PTR_ARGS(FixedRegionPlaneBoundaryCondition<2>, p_bc4, (&cell_population, point, normal));
		simulator.AddCellPopulationBoundaryCondition(p_bc4);

		simulator.Solve();
	}

};
Prev
Next