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

Analysis of fission in epithelial layer

Introduction

In this test we show how Chaste can be used to analyse the incidence of buckling in a layer of epithelial cells. Details of the computational model can be found in Almet 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 "DifferentiatedCellProliferativeType.hpp" //Non-proliferative cell type
#include "TransitCellProliferativeType.hpp" //Proliferative cell type
#include "WildtypeCellMutationState.hpp" //Standard wildtype state for all cells
#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 "HardCellMutationState.hpp" //Mutation class that defines hard 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 TestCryptFissionSweepsLiteratePaper : public AbstractCellBasedTestSuite
{
public:
	void TestEpithelialLayerSweeps() 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 = end_time/dt; //Set sampling timestep multiple (we only need the final state for analysis purposes

		//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 relative adhesiveness of hard cells. The values used in the paper were 1.0, 1.3 and 2.0
		double hard_cell_drag_multiplier = 1.0;

		//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 layer_radius = circle_radius + 2.0; //Radius of the layer of cells. This isn't the actual radius, just has to be large enough for later
		assert((layer_radius <= cells_across)&&(layer_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);
			}
		}

		//This is for the purpose of tracking cells. Weird computer science stuff.
		boost::shared_ptr<AbstractCellProperty> p_diff_type = CellPropertyRegistry::Instance()->Get<DifferentiatedCellProliferativeType>();
		boost::shared_ptr<AbstractCellProperty> p_soft_type = CellPropertyRegistry::Instance()->Get<TransitCellProliferativeType>();
		boost::shared_ptr<AbstractCellProperty> p_hard_state = CellPropertyRegistry::Instance()->Get<HardCellMutationState>();
		boost::shared_ptr<AbstractCellProperty> p_wildtype_state = CellPropertyRegistry::Instance()->Get<WildTypeCellMutationState>();

Reseed the random number generator via a command line argument

		unsigned index = CommandLineArguments::Instance()->GetUnsignedCorrespondingToOption("-myintval");
		RandomNumberGenerator::Instance()->Reseed(100*index);

Perform a parameter sweep over the stiffness ratio ond target proportions

		for (double hstiff = 0.0; hstiff < 9.0; hstiff ++) //Sweep over 3 <= \mu_H/\mu_S <= 5 in increments of 0.25
		{
			for (double tpropn = 0.0; tpropn < 15.0; tpropn ++) //Sweep over 0.25 <= p <= 0.95 in increments of 0.05
			{

				double stiffness_ratio = 3.0 + 0.25*hstiff;
				double target_proportion = 0.25 + 0.05*tpropn;

				//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_wildtype_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 layer 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(layer_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_soft_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 layer to narrow down our search

					if (pow(x-circle_centre[0],2) + pow(y-circle_centre[1],2) <= pow(layer_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_nonepithelial_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_nonepithelial_nodes = true;
											break; 				// This should break out of the inner for loop
										}
									}
								}
							}

							if(!element_contains_nonepithelial_nodes)
							{
								cell_iter->Kill();
							}
						}
					}
				}

Randomly assign cells in the layer to be hard 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 layer, we draw a random number and assign cells to be soft 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_hard_state);
					}
				}

Define the simulation class.

				OffLatticeSimulation<2> simulator(cell_population);

Set the drag constant of hard cells relative to soft cells

				double normal_damping_constant = cell_population.GetDampingConstantNormal();
				cell_population.SetDampingConstantMutant(hard_cell_drag_multiplier*normal_damping_constant);

				//Simulation pre-amble

				//Set output directory
				std::stringstream out;
				out << "Drag_" << hard_cell_drag_multiplier;
				out << "/BM_" << bm_force << "_TC_" << target_curvature;
				std::stringstream params;
				params << "SRatio_" << stiffness_ratio << "/TPropn_" << target_proportion;
				std::stringstream run;
				run << index;
				std::string output_directory = "CryptFissionSweepsLiteratePaper/" + out.str() + "/"
						+ run.str() + "/" + params.str();
				simulator.SetOutputDirectory(output_directory);

				simulator.SetDt(dt);
				simulator.SetSamplingTimestepMultiple(sampling_timestep); //Sample the simulation at every hour
				simulator.SetEndTime(end_time); //Hopefully this is long enough for a steady state

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->SetHardCellStiffnessRatio(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();

We now collect statistics on the positions and areas of the cells in the epithelial layer that are used for subsequent analysis.

We first need the epithelial cells ordered around the layer.

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

				c_vector<double,2> centre_of_mass; //Initialise centre of mass

				unsigned num_epithelial_cells = 0; //Initialise count

				for (AbstractCellPopulation<2>::Iterator cell_iter = simulator.rGetCellPopulation().Begin();
						cell_iter != simulator.rGetCellPopulation().End(); ++cell_iter)
				{
					//Only consider epithelial cells
					if(!cell_iter->GetCellProliferativeType()->IsType<DifferentiatedCellProliferativeType>() )
					{
						//Get node inex
						unsigned node_index = simulator.rGetCellPopulation().GetLocationIndexUsingCell(*cell_iter);

						//Get location of cell
						c_vector<double,2> cell_location = simulator.rGetCellPopulation().GetLocationOfCellCentre(*cell_iter);

						indices_of_cells_in_layer.push_back(node_index); //Add node index

						centre_of_mass += cell_location; //Update the centre of mass

						num_epithelial_cells += 1; //Update count of epithelial cells in layer
					}
				}

				//Normalise centre_of_mass
				centre_of_mass /= num_epithelial_cells;

				//Initialise vector of angles and the cell indices
				std::vector<std::pair<double, unsigned> > angles_and_cell_indices;

Get the angle of each epithelial cell with respect to the centre of mass. That way we can sort the list in counter-clockwise order before we create the file.

				for (unsigned k = 0; k < indices_of_cells_in_layer.size(); k++)
				{
					unsigned cell_index = indices_of_cells_in_layer[k]; //Get the node index

					CellPtr cell = simulator.rGetCellPopulation().GetCellUsingLocationIndex(cell_index); //Get the cell

					//Get the cell location
					double x = simulator.rGetCellPopulation().GetLocationOfCellCentre(cell)[0];
					double y = simulator.rGetCellPopulation().GetLocationOfCellCentre(cell)[1];

					//Make the pair of the angle and cell
					std::pair<double, unsigned> angle_cell;

					//Get point relative to the circle centre
					double rel_x = x - centre_of_mass[0];
					double rel_y = y - centre_of_mass[1];

					double circle_angle = atan(rel_y/rel_x); //Get initial angle argument

					if (rel_x<0.0) //If the point is in the second quadrant or third quadrant
					{
						circle_angle += M_PI;
					}
					else if ((rel_x>=0.0)&&(rel_y<0.0)) //Fourth quadrant
					{
						circle_angle += 2*M_PI;
					}

					angle_cell = std::make_pair(circle_angle, cell_index);

					angles_and_cell_indices.push_back(angle_cell); //Add the angle and node index
				}

				//Sort the vector of cell indices by their angle
				std::sort(angles_and_cell_indices.begin(), angles_and_cell_indices.end());

Create file to store the location of the cells in the layer and their cell type.

				OutputFileHandler results_handler(output_directory + "/", false); //Create output file handler

				std::string cells_in_layer_filename = "cellpositionsandareasinlayer.dat";
				out_stream cells_in_layer_file = results_handler.OpenOutputFile(cells_in_layer_filename);

Add each cell location and area to the file, along with a 1 if the cell is a hard cell and 0 otherwise.

				for (unsigned k=0; k<angles_and_cell_indices.size(); k++)
				{
					std::pair<double, unsigned> angle_and_cell_index = angles_and_cell_indices[k]; //Get pair
					unsigned cell_index = angle_and_cell_index.second; //Get cell index

					CellPtr cell = simulator.rGetCellPopulation().GetCellUsingLocationIndex(cell_index); //Get the cell

					//Get the area
					double cell_area = simulator.rGetCellPopulation().GetVolumeOfCell(cell);

					//Get the cell location
					double x = simulator.rGetCellPopulation().GetLocationOfCellCentre(cell)[0];
					double y = simulator.rGetCellPopulation().GetLocationOfCellCentre(cell)[1];

					unsigned is_hard_cell = 1;

					if(!cell->GetMutationState()->IsType<HardCellMutationState>() )
					{
						is_hard_cell = 0;
					}

					*cells_in_layer_file << x << "\t" << y << "\t" << cell_area << is_hard_cell << "\n";
				}

				//Tidying up
				SimulationTime::Instance()->Destroy();
				SimulationTime::Instance()->SetStartTime(0.0);
			}
		}
	}

};

Code

The full code is given below

File name TestCryptFissionSweepsLiteratePaper.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 "DifferentiatedCellProliferativeType.hpp" //Non-proliferative cell type
#include "TransitCellProliferativeType.hpp" //Proliferative cell type
#include "WildtypeCellMutationState.hpp" //Standard wildtype state for all cells
#include "FakePetscSetup.hpp" //Forbids tests running in parallel

#include "StochasticTargetProportionBasedCellCycleModel.hpp" //Asymmetric-division-based cell cycle model
#include "HardCellMutationState.hpp" //Mutation class that defines hard 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 TestCryptFissionSweepsLiteratePaper : public AbstractCellBasedTestSuite
{
public:
	void TestEpithelialLayerSweeps() throw(Exception)
	{
		//Simulation time parameters
		double dt = 0.005; //Set dt
		double end_time = 100.0; //Set end time
		double sampling_timestep = end_time/dt; //Set sampling timestep multiple (we only need the final state for analysis purposes

		//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 relative adhesiveness of hard cells. The values used in the paper were 1.0, 1.3 and 2.0
		double hard_cell_drag_multiplier = 1.0;

		//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 layer_radius = circle_radius + 2.0; //Radius of the layer of cells. This isn't the actual radius, just has to be large enough for later
		assert((layer_radius <= cells_across)&&(layer_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);
			}
		}

		//This is for the purpose of tracking cells. Weird computer science stuff.
		boost::shared_ptr<AbstractCellProperty> p_diff_type = CellPropertyRegistry::Instance()->Get<DifferentiatedCellProliferativeType>();
		boost::shared_ptr<AbstractCellProperty> p_soft_type = CellPropertyRegistry::Instance()->Get<TransitCellProliferativeType>();
		boost::shared_ptr<AbstractCellProperty> p_hard_state = CellPropertyRegistry::Instance()->Get<HardCellMutationState>();
		boost::shared_ptr<AbstractCellProperty> p_wildtype_state = CellPropertyRegistry::Instance()->Get<WildTypeCellMutationState>();

		unsigned index = CommandLineArguments::Instance()->GetUnsignedCorrespondingToOption("-myintval");
		RandomNumberGenerator::Instance()->Reseed(100*index);

		for (double hstiff = 0.0; hstiff < 9.0; hstiff ++) //Sweep over 3 <= \mu_H/\mu_S <= 5 in increments of 0.25
		{
			for (double tpropn = 0.0; tpropn < 15.0; tpropn ++) //Sweep over 0.25 <= p <= 0.95 in increments of 0.05
			{

				double stiffness_ratio = 3.0 + 0.25*hstiff;
				double target_proportion = 0.25 + 0.05*tpropn;

				//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_wildtype_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(layer_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_soft_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 layer to narrow down our search

					if (pow(x-circle_centre[0],2) + pow(y-circle_centre[1],2) <= pow(layer_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_nonepithelial_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_nonepithelial_nodes = true;
											break; 				// This should break out of the inner for loop
										}
									}
								}
							}

							if(!element_contains_nonepithelial_nodes)
							{
								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_hard_state);
					}
				}

				OffLatticeSimulation<2> simulator(cell_population);

				double normal_damping_constant = cell_population.GetDampingConstantNormal();
				cell_population.SetDampingConstantMutant(hard_cell_drag_multiplier*normal_damping_constant);

				//Simulation pre-amble

				//Set output directory
				std::stringstream out;
				out << "Drag_" << hard_cell_drag_multiplier;
				out << "/BM_" << bm_force << "_TC_" << target_curvature;
				std::stringstream params;
				params << "SRatio_" << stiffness_ratio << "/TPropn_" << target_proportion;
				std::stringstream run;
				run << index;
				std::string output_directory = "CryptFissionSweepsLiteratePaper/" + out.str() + "/"
						+ run.str() + "/" + params.str();
				simulator.SetOutputDirectory(output_directory);

				simulator.SetDt(dt);
				simulator.SetSamplingTimestepMultiple(sampling_timestep); //Sample the simulation at every hour
				simulator.SetEndTime(end_time); //Hopefully this is long enough for a steady state

				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->SetHardCellStiffnessRatio(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();

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

				c_vector<double,2> centre_of_mass; //Initialise centre of mass

				unsigned num_epithelial_cells = 0; //Initialise count

				for (AbstractCellPopulation<2>::Iterator cell_iter = simulator.rGetCellPopulation().Begin();
						cell_iter != simulator.rGetCellPopulation().End(); ++cell_iter)
				{
					//Only consider epithelial cells
					if(!cell_iter->GetCellProliferativeType()->IsType<DifferentiatedCellProliferativeType>() )
					{
						//Get node inex
						unsigned node_index = simulator.rGetCellPopulation().GetLocationIndexUsingCell(*cell_iter);

						//Get location of cell
						c_vector<double,2> cell_location = simulator.rGetCellPopulation().GetLocationOfCellCentre(*cell_iter);

						indices_of_cells_in_layer.push_back(node_index); //Add node index

						centre_of_mass += cell_location; //Update the centre of mass

						num_epithelial_cells += 1; //Update count of epithelial cells in layer
					}
				}

				//Normalise centre_of_mass
				centre_of_mass /= num_epithelial_cells;

				//Initialise vector of angles and the cell indices
				std::vector<std::pair<double, unsigned> > angles_and_cell_indices;

				for (unsigned k = 0; k < indices_of_cells_in_layer.size(); k++)
				{
					unsigned cell_index = indices_of_cells_in_layer[k]; //Get the node index

					CellPtr cell = simulator.rGetCellPopulation().GetCellUsingLocationIndex(cell_index); //Get the cell

					//Get the cell location
					double x = simulator.rGetCellPopulation().GetLocationOfCellCentre(cell)[0];
					double y = simulator.rGetCellPopulation().GetLocationOfCellCentre(cell)[1];

					//Make the pair of the angle and cell
					std::pair<double, unsigned> angle_cell;

					//Get point relative to the circle centre
					double rel_x = x - centre_of_mass[0];
					double rel_y = y - centre_of_mass[1];

					double circle_angle = atan(rel_y/rel_x); //Get initial angle argument

					if (rel_x<0.0) //If the point is in the second quadrant or third quadrant
					{
						circle_angle += M_PI;
					}
					else if ((rel_x>=0.0)&&(rel_y<0.0)) //Fourth quadrant
					{
						circle_angle += 2*M_PI;
					}

					angle_cell = std::make_pair(circle_angle, cell_index);

					angles_and_cell_indices.push_back(angle_cell); //Add the angle and node index
				}

				//Sort the vector of cell indices by their angle
				std::sort(angles_and_cell_indices.begin(), angles_and_cell_indices.end());

				OutputFileHandler results_handler(output_directory + "/", false); //Create output file handler

				std::string cells_in_layer_filename = "cellpositionsandareasinlayer.dat";
				out_stream cells_in_layer_file = results_handler.OpenOutputFile(cells_in_layer_filename);

				for (unsigned k=0; k<angles_and_cell_indices.size(); k++)
				{
					std::pair<double, unsigned> angle_and_cell_index = angles_and_cell_indices[k]; //Get pair
					unsigned cell_index = angle_and_cell_index.second; //Get cell index

					CellPtr cell = simulator.rGetCellPopulation().GetCellUsingLocationIndex(cell_index); //Get the cell

					//Get the area
					double cell_area = simulator.rGetCellPopulation().GetVolumeOfCell(cell);

					//Get the cell location
					double x = simulator.rGetCellPopulation().GetLocationOfCellCentre(cell)[0];
					double y = simulator.rGetCellPopulation().GetLocationOfCellCentre(cell)[1];

					unsigned is_hard_cell = 1;

					if(!cell->GetMutationState()->IsType<HardCellMutationState>() )
					{
						is_hard_cell = 0;
					}

					*cells_in_layer_file << x << "\t" << y << "\t" << cell_area << is_hard_cell << "\n";
				}

				//Tidying up
				SimulationTime::Instance()->Destroy();
				SimulationTime::Instance()->SetStartTime(0.0);
			}
		}
	}

};