From 609afbaa9eb642fa39b633fe900ae0423d877f7e Mon Sep 17 00:00:00 2001 From: Matteo Date: Tue, 7 Feb 2017 15:33:28 +0000 Subject: [PATCH] first commit --- CMakeLists.txt | 6 +- apps/CMakeLists.txt | 2 +- src/ConstantTargetAreaModifier.cpp | 73 ++++++++++++ src/ConstantTargetAreaModifier.hpp | 94 +++++++++++++++ src/MatteoCellCycleModel.cpp | 81 +++++++++++++ src/MatteoCellCycleModel.hpp | 178 +++++++++++++++++++++++++++++ src/MatteoForce.cpp | 64 +++++++++++ src/MatteoForce.hpp | 73 ++++++++++++ src/MatteoModifier.cpp | 138 ++++++++++++++++++++++ src/MatteoModifier.hpp | 118 +++++++++++++++++++ src/MatteoSrnModel.cpp | 137 ++++++++++++++++++++++ src/MatteoSrnModel.hpp | 148 ++++++++++++++++++++++++ test/CMakeLists.txt | 2 +- test/ContinuousTestPack.txt | 4 +- test/TestOptogenetics.hpp | 100 ++++++++++++++++ test/Testmatteo.hpp | 114 ++++++++++++++++++ 16 files changed, 1326 insertions(+), 6 deletions(-) create mode 100644 src/ConstantTargetAreaModifier.cpp create mode 100644 src/ConstantTargetAreaModifier.hpp create mode 100644 src/MatteoCellCycleModel.cpp create mode 100644 src/MatteoCellCycleModel.hpp create mode 100644 src/MatteoForce.cpp create mode 100644 src/MatteoForce.hpp create mode 100644 src/MatteoModifier.cpp create mode 100644 src/MatteoModifier.hpp create mode 100644 src/MatteoSrnModel.cpp create mode 100644 src/MatteoSrnModel.hpp create mode 100644 test/TestOptogenetics.hpp create mode 100644 test/Testmatteo.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 37e8007..7b53ad2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,7 +3,7 @@ # You can set which Chaste components (or other projects) your project depends on by editing the # find_package() call for Chaste. E.g. -# find_package(Chaste COMPONENTS cell_based) + find_package(Chaste COMPONENTS cell_based) # for a project just using the cell_based component (and its dependencies), or # find_package(Chaste COMPONENTS heart lung) # for a project combining heart & lung simulations. @@ -11,7 +11,7 @@ # Note that the order in which components are specified does not matter. # Here we just depend on core components (nothing application-specific). -find_package(Chaste COMPONENTS continuum_mechanics global io linalg mesh ode pde) +#find_package(Chaste COMPONENTS cell_based global io linalg mesh ode pde) # Alternatively, to specify a Chaste installation directory use a line like that below. # This is needed if your project is not contained in the projects folder within a Chaste source tree. @@ -19,4 +19,4 @@ find_package(Chaste COMPONENTS continuum_mechanics global io linalg mesh ode pde # Change the project name in the line below to match the folder this file is in, # i.e. the name of your project. -chaste_do_project(template_project) +chaste_do_project(matteo_chaste) diff --git a/apps/CMakeLists.txt b/apps/CMakeLists.txt index 439d427..0034b25 100644 --- a/apps/CMakeLists.txt +++ b/apps/CMakeLists.txt @@ -3,4 +3,4 @@ # This allows us to build executables (apps) using the project. # Simply set the project name in the line below. -chaste_do_apps_project(template_project) +chaste_do_apps_project(matteo_chaste) diff --git a/src/ConstantTargetAreaModifier.cpp b/src/ConstantTargetAreaModifier.cpp new file mode 100644 index 0000000..93694e5 --- /dev/null +++ b/src/ConstantTargetAreaModifier.cpp @@ -0,0 +1,73 @@ +/* + +Copyright (c) 2005-2016, University of Oxford. +All rights reserved. + +University of Oxford means the Chancellor, Masters and Scholars of the +University of Oxford, having an administrative office at Wellington +Square, Oxford OX1 2JD, UK. + +This file is part of Chaste. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + * Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + * Neither the name of the University of Oxford nor the names of its + contributors may be used to endorse or promote products derived from this + software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE +GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT +OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +*/ + +#include "ConstantTargetAreaModifier.hpp" + +template +ConstantTargetAreaModifier::ConstantTargetAreaModifier() + : AbstractTargetAreaModifier() +{ +} + +template +ConstantTargetAreaModifier::~ConstantTargetAreaModifier() +{ +} + +template +void ConstantTargetAreaModifier::UpdateTargetAreaOfCell(CellPtr pCell) +{ + // Get target area A of a healthy cell in S, G2 or M phase + double cell_target_area = this->mReferenceTargetArea; + + // Set cell data + pCell->GetCellData()->SetItem("target area", cell_target_area); +} + +template +void ConstantTargetAreaModifier::OutputSimulationModifierParameters(out_stream& rParamsFile) +{ + // No parameters to output, so just call method on direct parent class + AbstractTargetAreaModifier::OutputSimulationModifierParameters(rParamsFile); +} + +// Explicit instantiation +template class ConstantTargetAreaModifier<1>; +template class ConstantTargetAreaModifier<2>; +template class ConstantTargetAreaModifier<3>; + +// Serialization for Boost >= 1.36 +#include "SerializationExportWrapperForCpp.hpp" +EXPORT_TEMPLATE_CLASS_SAME_DIMS(ConstantTargetAreaModifier) diff --git a/src/ConstantTargetAreaModifier.hpp b/src/ConstantTargetAreaModifier.hpp new file mode 100644 index 0000000..f53ae37 --- /dev/null +++ b/src/ConstantTargetAreaModifier.hpp @@ -0,0 +1,94 @@ +/* + +Copyright (c) 2005-2016, University of Oxford. +All rights reserved. + +University of Oxford means the Chancellor, Masters and Scholars of the +University of Oxford, having an administrative office at Wellington +Square, Oxford OX1 2JD, UK. + +This file is part of Chaste. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + * Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + * Neither the name of the University of Oxford nor the names of its + contributors may be used to endorse or promote products derived from this + software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE +GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT +OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +*/ + +#ifndef CONSTANTTARGETAREAMODIFIER_HPP_ +#define CONSTANTTARGETAREAMODIFIER_HPP_ + +#include "ChasteSerialization.hpp" +#include "AbstractTargetAreaModifier.hpp" + +/** + * A modifier class in which the target area property of each cell is held constant over time. + */ +template +class ConstantTargetAreaModifier : public AbstractTargetAreaModifier +{ + /** Needed for serialization. */ + friend class boost::serialization::access; + /** + * Boost Serialization method for archiving/checkpointing. + * Archives the object and its member variables. + * + * @param archive The boost archive. + * @param version The current version of this class. + */ + template + void serialize(Archive & archive, const unsigned int version) + { + archive & boost::serialization::base_object >(*this); + } + +public: + + /** + * Default constructor. + */ + ConstantTargetAreaModifier(); + + /** + * Destructor. + */ + virtual ~ConstantTargetAreaModifier(); + + /** + * Overridden UpdateTargetAreaOfCell() method. + * + * @param pCell pointer to the cell + */ + virtual void UpdateTargetAreaOfCell(const CellPtr pCell); + + /** + * Overridden OutputSimulationModifierParameters() method. + * Output any simulation modifier parameters to file. + * + * @param rParamsFile the file stream to which the parameters are output + */ + void OutputSimulationModifierParameters(out_stream& rParamsFile); +}; + +#include "SerializationExportWrapper.hpp" +EXPORT_TEMPLATE_CLASS_SAME_DIMS(ConstantTargetAreaModifier) + +#endif /*CONSTANTTARGETAREAMODIFIER_HPP_*/ diff --git a/src/MatteoCellCycleModel.cpp b/src/MatteoCellCycleModel.cpp new file mode 100644 index 0000000..ea7885c --- /dev/null +++ b/src/MatteoCellCycleModel.cpp @@ -0,0 +1,81 @@ + +#include "MatteoCellCycleModel.hpp" +#include "DifferentiatedCellProliferativeType.hpp" +#include "Debug.hpp" + +MatteoCellCycleModel::MatteoCellCycleModel() + : AbstractSimpleCellCycleModel(), + mMinCellCycleDuration(12.0), // Hours + mMaxCellCycleDuration(14.0) // Hours +{ +} + +/*called every time step and tells if a cell is aready to divide (using the cell cycle model)*/ + +bool MatteoCellCycleModel::ReadyToDivide() +{ + bool ready = mpCell->GetCellData()->GetItem("divide"); + mpCell->GetCellData()->SetItem("divide", 0); + return ready; + + +} + +MatteoCellCycleModel::MatteoCellCycleModel(const MatteoCellCycleModel& rModel) + : AbstractSimpleCellCycleModel(rModel), + mMinCellCycleDuration(rModel.mMinCellCycleDuration), + mMaxCellCycleDuration(rModel.mMaxCellCycleDuration) +{ +} + +AbstractCellCycleModel* MatteoCellCycleModel::CreateCellCycleModel() +{ + return new MatteoCellCycleModel(*this); +} + +void MatteoCellCycleModel::SetCellCycleDuration() +{ +} + +double MatteoCellCycleModel::GetMinCellCycleDuration() +{ + return mMinCellCycleDuration; +} + +void MatteoCellCycleModel::SetMinCellCycleDuration(double minCellCycleDuration) +{ + mMinCellCycleDuration = minCellCycleDuration; +} + +double MatteoCellCycleModel::GetMaxCellCycleDuration() +{ + return mMaxCellCycleDuration; +} + +void MatteoCellCycleModel::SetMaxCellCycleDuration(double maxCellCycleDuration) +{ + mMaxCellCycleDuration = maxCellCycleDuration; +} + +double MatteoCellCycleModel::GetAverageTransitCellCycleTime() +{ + return 0.5*(mMinCellCycleDuration + mMaxCellCycleDuration); +} + +double MatteoCellCycleModel::GetAverageStemCellCycleTime() +{ + return 0.5*(mMinCellCycleDuration + mMaxCellCycleDuration); +} + +void MatteoCellCycleModel::OutputCellCycleModelParameters(out_stream& rParamsFile) +{ + *rParamsFile << "\t\t\t" << mMinCellCycleDuration << "\n"; + *rParamsFile << "\t\t\t" << mMaxCellCycleDuration << "\n"; + + // Call method on direct parent class + AbstractSimpleCellCycleModel::OutputCellCycleModelParameters(rParamsFile); +} + +// Serialization for Boost >= 1.36 +#include "SerializationExportWrapperForCpp.hpp" +CHASTE_CLASS_EXPORT(MatteoCellCycleModel) diff --git a/src/MatteoCellCycleModel.hpp b/src/MatteoCellCycleModel.hpp new file mode 100644 index 0000000..a6ff8fd --- /dev/null +++ b/src/MatteoCellCycleModel.hpp @@ -0,0 +1,178 @@ +/* + +Copyright (c) 2005-2016, University of Oxford. +All rights reserved. + +University of Oxford means the Chancellor, Masters and Scholars of the +University of Oxford, having an administrative office at Wellington +Square, Oxford OX1 2JD, UK. + +This file is part of Chaste. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + * Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + * Neither the name of the University of Oxford nor the names of its + contributors may be used to endorse or promote products derived from this + software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE +GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT +OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +*/ + +#ifndef MatteoCellCycleModel_HPP_ +#define MatteoCellCycleModel_HPP_ + +#include "AbstractSimpleCellCycleModel.hpp" +#include "RandomNumberGenerator.hpp" + +/** + * A stochastic cell-cycle model where cells divide with a stochastic cell cycle duration + * with the length of the cell cycle drawn from a uniform distribution + * on [mMinCellCycleDuration, mMaxCellCycleDuration]. + * + * If the cell is differentiated, then the cell cycle duration is set to be infinite, + * so that the cell will never divide. + */ +class MatteoCellCycleModel : public AbstractSimpleCellCycleModel +{ + friend class TestSimpleCellCycleModels; + +private: + + /** + * The minimum cell cycle duration. Used to define the uniform distribution. + * Defaults to 12 hours. + */ + double mMinCellCycleDuration; + + /** + * The maximum cell cycle duration. Used to define the uniform distribution. + * Defaults to 14 hours. + */ + double mMaxCellCycleDuration; + + /** Needed for serialization. */ + friend class boost::serialization::access; + /** + * Archive the cell-cycle model and random number generator, never used directly - boost uses this. + * + * @param archive the archive + * @param version the current version of this class + */ + template + void serialize(Archive & archive, const unsigned int version) + { + archive & boost::serialization::base_object(*this); + + // Make sure the RandomNumberGenerator singleton gets saved too + SerializableSingleton* p_wrapper = RandomNumberGenerator::Instance()->GetSerializationWrapper(); + archive & p_wrapper; + archive & mMinCellCycleDuration; + archive & mMaxCellCycleDuration; + } + +protected: + + /** + * Protected copy-constructor for use by CreateCellCycleModel(). + * + * The only way for external code to create a copy of a cell cycle model + * is by calling that method, to ensure that a model of the correct subclass is created. + * This copy-constructor helps subclasses to ensure that all member variables are correctly copied when this happens. + * + * This method is called by child classes to set member variables for a daughter cell upon cell division. + * Note that the parent cell cycle model will have had ResetForDivision() called just before CreateCellCycleModel() is called, + * so performing an exact copy of the parent is suitable behaviour. Any daughter-cell-specific initialisation + * can be done in InitialiseDaughterCell(). + * + * @param rModel the cell cycle model to copy. + */ + MatteoCellCycleModel(const MatteoCellCycleModel& rModel); + +public: + + /** + * Constructor - just a default, mBirthTime is set in the AbstractCellCycleModel class. + * mG1Duration is set very high, it is set for the individual cells when InitialiseDaughterCell is called + */ + MatteoCellCycleModel(); + + /** + * Overridden SetCellCycleDuration Method to add stochastic cell cycle times + */ + void SetCellCycleDuration(); + + bool ReadyToDivide(); + /** + * Overridden builder method to create new copies of + * this cell-cycle model. + * + * @return new cell-cycle model + */ + AbstractCellCycleModel* CreateCellCycleModel(); + + /** + * @return mMinCellCycleDuration + */ + double GetMinCellCycleDuration(); + + /** + * Set mMinCellCycleDuration. + * + * @param minCellCycleDuration + */ + void SetMinCellCycleDuration(double minCellCycleDuration); + + /** + * @return mMaxCellCycleDuration + */ + double GetMaxCellCycleDuration(); + + /** + * Set mMaxCellCycleDuration. + * + * @param maxCellCycleDuration + */ + void SetMaxCellCycleDuration(double maxCellCycleDuration); + + /** + * Overridden GetAverageTransitCellCycleTime() method. + * + * @return the average of mMinCellCycleDuration and mMaxCellCycleDuration + */ + double GetAverageTransitCellCycleTime(); + + /** + * Overridden GetAverageStemCellCycleTime() method. + * + * @return the average of mMinCellCycleDuration and mMaxCellCycleDuration + */ + double GetAverageStemCellCycleTime(); + + /** + * Overridden OutputCellCycleModelParameters() method. + * + * @param rParamsFile the file stream to which the parameters are output + */ + virtual void OutputCellCycleModelParameters(out_stream& rParamsFile); +}; + +#include "SerializationExportWrapper.hpp" +// Declare identifier for the serializer +CHASTE_CLASS_EXPORT(MatteoCellCycleModel) + +#endif /*MatteoCellCycleModel_HPP_*/ diff --git a/src/MatteoForce.cpp b/src/MatteoForce.cpp new file mode 100644 index 0000000..590ca0d --- /dev/null +++ b/src/MatteoForce.cpp @@ -0,0 +1,64 @@ + + +#include "MatteoForce.hpp" + +template +MatteoForce::MatteoForce() + : FarhadifarForce() + { +} + +template +MatteoForce::~MatteoForce() +{ +} + + +template +double MatteoForce::GetLineTensionParameter(Node* pNodeA, Node* pNodeB, VertexBasedCellPopulation& rVertexCellPopulation) +{ + // Find the indices of the elements owned by each node + std::set elements_containing_nodeA = pNodeA->rGetContainingElementIndices(); + std::set elements_containing_nodeB = pNodeB->rGetContainingElementIndices(); + + // Find common elements + std::set shared_elements; + std::set_intersection(elements_containing_nodeA.begin(), + elements_containing_nodeA.end(), + elements_containing_nodeB.begin(), + elements_containing_nodeB.end(), + std::inserter(shared_elements, shared_elements.begin())); + + // Check that the nodes have a common edge + assert(!shared_elements.empty()); + + // Since each internal edge is visited twice in the loop above, we have to use half the line tension parameter + // for each visit. + //double line_tension_parameter_in_calculation = this->GetLineTensionParameter()/2.0; + double line_tension_parameter_in_calculation = 0.12; + // If the edge corresponds to a single element, then the cell is on the boundary + if (shared_elements.size() == 1) + { + // line_tension_parameter_in_calculation = this->GetBoundaryLineTensionParameter(); + line_tension_parameter_in_calculation = 0.12; + } + + return line_tension_parameter_in_calculation; +} + + +template +void MatteoForce::OutputForceParameters(out_stream& rParamsFile) +{ + // Call method on direct parent class + FarhadifarForce::OutputForceParameters(rParamsFile); +} + +// Explicit instantiation +template class MatteoForce<1>; +template class MatteoForce<2>; +template class MatteoForce<3>; + +// Serialization for Boost >= 1.36 +#include "SerializationExportWrapperForCpp.hpp" +EXPORT_TEMPLATE_CLASS_SAME_DIMS(MatteoForce) diff --git a/src/MatteoForce.hpp b/src/MatteoForce.hpp new file mode 100644 index 0000000..6317af2 --- /dev/null +++ b/src/MatteoForce.hpp @@ -0,0 +1,73 @@ + +#ifndef MatteoFORCE_HPP_ +#define MatteoFORCE_HPP_ + +#include "ChasteSerialization.hpp" +#include +#include "Exception.hpp" + +#include "FarhadifarForce.hpp" +#include "VertexBasedCellPopulation.hpp" + +#include + +/** + * A force class for use in Vertex-based simulations. This force is based on the + * Energy function proposed by Matteo et al in Curr. Biol., 2007, 17, 2095-2104. + */ + + +template +class MatteoForce : public FarhadifarForce +{ +friend class TestForces; + +private: + + friend class boost::serialization::access; + /** + * Boost Serialization method for archiving/checkpointing. + * Archives the object and its member variables. + * + * @param archive The boost archive. + * @param version The current version of this class. + */ + template + void serialize(Archive & archive, const unsigned int version) + { + archive & boost::serialization::base_object >(*this); + + } + + +public: + + /** + * Constructor. + */ + MatteoForce(); + + /** + * Destructor. + */ + virtual ~MatteoForce(); + + + /** + * Get the line tension parameter for the edge between two given nodes. + * + * @param pNodeA one node + * @param pNodeB the other node + * @param rVertexCellPopulation reference to the cell population + * + * @return the line tension parameter for this edge. + */ + virtual double GetLineTensionParameter(Node* pNodeA, Node* pNodeB, VertexBasedCellPopulation& rVertexCellPopulation); + + void OutputForceParameters(out_stream& rParamsFile); +}; + +#include "SerializationExportWrapper.hpp" +EXPORT_TEMPLATE_CLASS_SAME_DIMS(MatteoForce) + +#endif /*MatteoFORCE_HPP_*/ diff --git a/src/MatteoModifier.cpp b/src/MatteoModifier.cpp new file mode 100644 index 0000000..a541a2a --- /dev/null +++ b/src/MatteoModifier.cpp @@ -0,0 +1,138 @@ + +#include "MatteoModifier.hpp" +#include "RandomNumberGenerator.hpp" +#include "CellLabel.hpp" +#include "Debug.hpp" + +template +MatteoModifier::MatteoModifier() + : AbstractCellBasedSimulationModifier() +{ +} + +template +MatteoModifier::~MatteoModifier() +{ +} + +template +void MatteoModifier::UpdateAtEndOfTimeStep(AbstractCellPopulation& rCellPopulation) +{ + + + //do every 10 units of time + if (fmod(SimulationTime::Instance()->GetTime(),10) == 0) + { + // Store each cell's current fitness + UpdateCellData(rCellPopulation); + + // Randomly pick one cell to divide + std::map map; + double prev = rCellPopulation.Begin()->GetCellData()->GetItem("fitness"); + for (typename AbstractCellPopulation::Iterator cell_iter = rCellPopulation.Begin(); + cell_iter != rCellPopulation.End(); + ++cell_iter) + { + if (*cell_iter == *(rCellPopulation.Begin())) + { + map[*cell_iter] = prev; + } + else + { + map[*cell_iter] = prev + cell_iter->GetCellData()->GetItem("fitness"); + prev +=cell_iter->GetCellData()->GetItem("fitness"); + } + } + double total_fitness = prev; + for (std::map::iterator iter = map.begin(); iter!=map.end(); ++iter) + { + iter->second /= total_fitness; + } + + double r = RandomNumberGenerator::Instance()->ranf(); + for (std::map::iterator iter = map.begin(); iter!=map.end(); ++iter) + { + std::map::iterator other_iter = iter; + ++other_iter; + if ((iter->second < r) && (other_iter->second >= r)) + { + TellCellToDivide(iter->first); + break; + } + } + } + + +} + +template +void MatteoModifier::TellCellToDivide(CellPtr pCell) +{ +pCell->GetCellData()->SetItem("divide", 1); +} + +template +void MatteoModifier::SetupSolve(AbstractCellPopulation& rCellPopulation, std::string outputDirectory) +{ +} + +template +void MatteoModifier::UpdateCellData(AbstractCellPopulation& rCellPopulation) +{ +double b = 10; +double c = 5; +double delta = 0.01; + // Iterate over cell population + for (typename AbstractCellPopulation::Iterator cell_iter = rCellPopulation.Begin(); + cell_iter != rCellPopulation.End(); + ++cell_iter) + { + bool is_defector = cell_iter->template HasCellProperty(); + double payoff = 0; + std::set neighbours = rCellPopulation.GetNeighbouringLocationIndices(*cell_iter); + for (std::set::iterator iter = neighbours.begin(); + iter != neighbours.end(); + ++iter) + { + CellPtr p_neighbour = rCellPopulation.GetCellUsingLocationIndex(*iter); + bool is_neighbour_defector = p_neighbour->template HasCellProperty(); + + double contribution = b - c; + if (is_defector && !is_neighbour_defector) + { + contribution = b; + } + else if (!is_defector && is_neighbour_defector) + { + contribution = -c; + } + else if (is_defector && is_neighbour_defector) + { + contribution = 0; + } + payoff += contribution; + } + // Get the fitness of this cell + double cell_fitness = pow(1 + delta, payoff); + + // Store the cell's fitness in CellData + cell_iter->GetCellData()->SetItem("fitness", cell_fitness); + } +} + +template +void MatteoModifier::OutputSimulationModifierParameters(out_stream& rParamsFile) +{ + // No parameters to output, so just call method on direct parent class + AbstractCellBasedSimulationModifier::OutputSimulationModifierParameters(rParamsFile); +} + +// Explicit instantiation +template class MatteoModifier<1>; +template class MatteoModifier<2>; +template class MatteoModifier<3>; + +// Serialization for Boost >= 1.36 +#include "SerializationExportWrapperForCpp.hpp" +EXPORT_TEMPLATE_CLASS_SAME_DIMS(MatteoModifier) + diff --git a/src/MatteoModifier.hpp b/src/MatteoModifier.hpp new file mode 100644 index 0000000..c154f09 --- /dev/null +++ b/src/MatteoModifier.hpp @@ -0,0 +1,118 @@ +/* + +Copyright (c) 2005-2016, University of Oxford. +All rights reserved. + +University of Oxford means the Chancellor, Masters and Scholars of the +University of Oxford, having an administrative office at Wellington +Square, Oxford OX1 2JD, UK. + +This file is part of Chaste. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + * Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + * Neither the name of the University of Oxford nor the names of its + contributors may be used to endorse or promote products derived from this + software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE +GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT +OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +*/ + +#ifndef MatteoModifier_HPP_ +#define MatteoModifier_HPP_ + +#include "ChasteSerialization.hpp" +#include + +#include "AbstractCellBasedSimulationModifier.hpp" + +/** + * A modifier class which at each simulation time step calculates the volume of each cell + * and stores it in in the CellData property as "volume". To be used in conjunction with + * contact inhibition cell cycle models. + */ +template +class MatteoModifier : public AbstractCellBasedSimulationModifier +{ + /** Needed for serialization. */ + friend class boost::serialization::access; + /** + * Boost Serialization method for archiving/checkpointing. + * Archives the object and its member variables. + * + * @param archive The boost archive. + * @param version The current version of this class. + */ + template + void serialize(Archive & archive, const unsigned int version) + { + archive & boost::serialization::base_object >(*this); + } + +public: + + /** + * Default constructor. + */ + MatteoModifier(); + + /** + * Destructor. + */ + virtual ~MatteoModifier(); + + /** + * Overridden UpdateAtEndOfTimeStep() method. + * + * Specify what to do in the simulation at the end of each time step. + * + * @param rCellPopulation reference to the cell population + */ + virtual void UpdateAtEndOfTimeStep(AbstractCellPopulation& rCellPopulation); + + /** + * Overridden SetupSolve() method. + * + * Specify what to do in the simulation before the start of the time loop. + * + * @param rCellPopulation reference to the cell population + * @param outputDirectory the output directory, relative to where Chaste output is stored + */ + virtual void SetupSolve(AbstractCellPopulation& rCellPopulation, std::string outputDirectory); + + /** + * Helper method to compute the volume of each cell in the population and store these in the CellData. + * + * @param rCellPopulation reference to the cell population + */ + void UpdateCellData(AbstractCellPopulation& rCellPopulation); + + /** + * Overridden OutputSimulationModifierParameters() method. + * Output any simulation modifier parameters to file. + * + * @param rParamsFile the file stream to which the parameters are output + */ + void OutputSimulationModifierParameters(out_stream& rParamsFile); +void TellCellToDivide(CellPtr pCell); +}; + +#include "SerializationExportWrapper.hpp" +EXPORT_TEMPLATE_CLASS_SAME_DIMS(MatteoModifier) + +#endif /*MatteoModifier_HPP_*/ diff --git a/src/MatteoSrnModel.cpp b/src/MatteoSrnModel.cpp new file mode 100644 index 0000000..ecba003 --- /dev/null +++ b/src/MatteoSrnModel.cpp @@ -0,0 +1,137 @@ +/* + +Copyright (c) 2005-2017, University of Oxford. +All rights reserved. + +University of Oxford means the Chancellor, Masters and Scholars of the +University of Oxford, having an administrative office at Wellington +Square, Oxford OX1 2JD, UK. + +This file is part of Chaste. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + * Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + * Neither the name of the University of Oxford nor the names of its + contributors may be used to endorse or promote products derived from this + software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE +GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT +OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +*/ + +#include "MatteoSrnModel.hpp" + +MatteoSrnModel::MatteoSrnModel(boost::shared_ptr pOdeSolver) + : AbstractOdeSrnModel(2, pOdeSolver) +{ + if (mpOdeSolver == boost::shared_ptr()) + { +#ifdef CHASTE_CVODE + mpOdeSolver = CellCycleModelOdeSolver::Instance(); + mpOdeSolver->Initialise(); + mpOdeSolver->SetMaxSteps(10000); +#else + mpOdeSolver = CellCycleModelOdeSolver::Instance(); + mpOdeSolver->Initialise(); + SetDt(0.001); +#endif //CHASTE_CVODE + } + assert(mpOdeSolver->IsSetUp()); +} + +MatteoSrnModel::MatteoSrnModel(const MatteoSrnModel& rModel) + : AbstractOdeSrnModel(rModel) +{ + /* + * Set each member variable of the new SRN model that inherits + * its value from the parent. + * + * Note 1: some of the new SRN model's member variables + * will already have been correctly initialized in its constructor. + * + * Note 2: one or more of the new SRN model's member variables + * may be set/overwritten as soon as InitialiseDaughterCell() is called on + * the new SRN model. + * + * Note 3: Only set the variables defined in this class. Variables defined + * in parent classes will be defined there. + */ + + assert(rModel.GetOdeSystem()); + SetOdeSystem(new DeltaNotchOdeSystem(rModel.GetOdeSystem()->rGetStateVariables())); +} + +AbstractSrnModel* MatteoSrnModel::CreateSrnModel() +{ + return new MatteoSrnModel(*this); +} + +void MatteoSrnModel::SimulateToCurrentTime() +{ + // Custom behaviour + UpdateMatteo(); + + // Run the ODE simulation as needed + AbstractOdeSrnModel::SimulateToCurrentTime(); +} + +void MatteoSrnModel::Initialise() +{ + AbstractOdeSrnModel::Initialise(new DeltaNotchOdeSystem); +} + +void MatteoSrnModel::UpdateMatteo() +{ + assert(mpOdeSystem != NULL); + assert(mpCell != NULL); + + double mean_delta = mpCell->GetCellData()->GetItem("mean delta"); + mpOdeSystem->SetParameter("Mean Delta", mean_delta); +} + +double MatteoSrnModel::GetNotch() +{ + assert(mpOdeSystem != NULL); + double notch = mpOdeSystem->rGetStateVariables()[0]; + return notch; +} + +double MatteoSrnModel::GetDelta() +{ + assert(mpOdeSystem != NULL); + double delta = mpOdeSystem->rGetStateVariables()[1]; + return delta; +} + +double MatteoSrnModel::GetMeanNeighbouringDelta() +{ + assert(mpOdeSystem != NULL); + double mean_neighbouring_delta = mpOdeSystem->GetParameter("Mean Delta"); + return mean_neighbouring_delta; +} + +void MatteoSrnModel::OutputSrnModelParameters(out_stream& rParamsFile) +{ + // No new parameters to output, so just call method on direct parent class + AbstractOdeSrnModel::OutputSrnModelParameters(rParamsFile); +} + +// Declare identifier for the serializer +#include "SerializationExportWrapperForCpp.hpp" +CHASTE_CLASS_EXPORT(MatteoSrnModel) +#include "CellCycleModelOdeSolverExportWrapper.hpp" +EXPORT_CELL_CYCLE_MODEL_ODE_SOLVER(MatteoSrnModel) diff --git a/src/MatteoSrnModel.hpp b/src/MatteoSrnModel.hpp new file mode 100644 index 0000000..40d5d51 --- /dev/null +++ b/src/MatteoSrnModel.hpp @@ -0,0 +1,148 @@ +/* + +Copyright (c) 2005-2017, University of Oxford. +All rights reserved. + +University of Oxford means the Chancellor, Masters and Scholars of the +University of Oxford, having an administrative office at Wellington +Square, Oxford OX1 2JD, UK. + +This file is part of Chaste. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + * Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + * Neither the name of the University of Oxford nor the names of its + contributors may be used to endorse or promote products derived from this + software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE +GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) +HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT +OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +*/ + +#ifndef MatteoSRNMODEL_HPP_ +#define MatteoSRNMODEL_HPP_ + +#include "ChasteSerialization.hpp" +#include + +#include "DeltaNotchOdeSystem.hpp" +#include "AbstractOdeSrnModel.hpp" + +/** + * A subclass of AbstractOdeSrnModel that includes a Delta-Notch ODE system in the sub-cellular reaction network. + * + * \todo #2752 document this class more thoroughly here + */ +class MatteoSrnModel : public AbstractOdeSrnModel +{ +private: + + /** Needed for serialization. */ + friend class boost::serialization::access; + /** + * Archive the cell-cycle model and member variables. + * + * @param archive the archive + * @param version the current version of this class + */ + template + void serialize(Archive & archive, const unsigned int version) + { + archive & boost::serialization::base_object(*this); + } + +protected: + /** + * Protected copy-constructor for use by CreateSrnModel. The only way for external code to create a copy of a SRN model + * is by calling that method, to ensure that a model of the correct subclass is created. + * This copy-constructor helps subclasses to ensure that all member variables are correctly copied when this happens. + * + * This method is called by child classes to set member variables for a daughter cell upon cell division. + * Note that the parent SRN model will have had ResetForDivision() called just before CreateSrnModel() is called, + * so performing an exact copy of the parent is suitable behaviour. Any daughter-cell-specific initialisation + * can be done in InitialiseDaughterCell(). + * + * @param rModel the SRN model to copy. + */ + MatteoSrnModel(const MatteoSrnModel& rModel); + +public: + + /** + * Default constructor calls base class. + * + * @param pOdeSolver An optional pointer to a cell-cycle model ODE solver object (allows the use of different ODE solvers) + */ + MatteoSrnModel(boost::shared_ptr pOdeSolver = boost::shared_ptr()); + + /** + * Overridden builder method to create new copies of + * this SRN model. + * + * @return a copy of the current SRN model. + */ + AbstractSrnModel* CreateSrnModel(); + + /** + * Initialise the cell-cycle model at the start of a simulation. + * + * This overridden method sets up a new Delta-Notch ODE system. + */ + void Initialise(); // override + + /** + * Overridden SimulateToTime() method for custom behaviour. + * + * \todo #2752 say what it does in this class + */ + void SimulateToCurrentTime(); + + /** + * Update the current levels of Delta and Notch in the cell. + */ + void UpdateMatteo(); + + /** + * @return the current Notch level in this cell. + */ + double GetNotch(); + + /** + * @return the current Delta level in this cell. + */ + double GetDelta(); + + /** + * @return the current level of Delta neighbouring the cell. + */ + double GetMeanNeighbouringDelta(); + + /** + * Output cell-cycle model parameters to file. + * + * @param rParamsFile the file stream to which the parameters are output + */ + void OutputSrnModelParameters(out_stream& rParamsFile); +}; + +// Declare identifier for the serializer +#include "SerializationExportWrapper.hpp" +CHASTE_CLASS_EXPORT(MatteoSrnModel) +#include "CellCycleModelOdeSolverExportWrapper.hpp" +EXPORT_CELL_CYCLE_MODEL_ODE_SOLVER(MatteoSrnModel) + +#endif /* MatteoSRNMODEL_HPP_ */ diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index c6bfa5d..fc3dff3 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -3,4 +3,4 @@ # This allows us to build tests of the project. # Simply set the project name in the line below. -chaste_do_test_project(template_project) +chaste_do_test_project(matteo_chaste) diff --git a/test/ContinuousTestPack.txt b/test/ContinuousTestPack.txt index 0e80dc8..7faa430 100644 --- a/test/ContinuousTestPack.txt +++ b/test/ContinuousTestPack.txt @@ -1 +1,3 @@ -TestHello.hpp \ No newline at end of file +TestHello.hpp +Testmatteo.hpp +TestOptogenetics.hpp diff --git a/test/TestOptogenetics.hpp b/test/TestOptogenetics.hpp new file mode 100644 index 0000000..c7e0ee6 --- /dev/null +++ b/test/TestOptogenetics.hpp @@ -0,0 +1,100 @@ + + +#ifndef TESTOPTOGENETICS_HPP_ +#define TESTOPTOGENETICS_HPP_ + + +#include +#include "CheckpointArchiveTypes.hpp" +#include "AbstractCellBasedTestSuite.hpp" +#include "HoneycombMeshGenerator.hpp" +#include "HoneycombVertexMeshGenerator.hpp" +#include "NodeBasedCellPopulation.hpp" +#include "OffLatticeSimulation.hpp" +#include "VertexBasedCellPopulation.hpp" +#include "MatteoForce.hpp" +#include "ConstantTargetAreaModifier.hpp" +#include "GeneralisedLinearSpringForce.hpp" +#include "DifferentiatedCellProliferativeType.hpp" +#include "WildTypeCellMutationState.hpp" +#include "SmartPointers.hpp" +#include "PetscSetupAndFinalize.hpp" +#include "NoCellCycleModel.hpp" + +/* + * The next header file defines a simple subcellular reaction network model that includes the functionality + * for solving each cell's Delta/Notch signalling ODE system at each time step, using information about neighbouring + * cells through the {{{CellData}}} class. + */ +#include "MatteoSrnModel.hpp" +/* + * The next header defines the simulation class modifier corresponding to the Delta-Notch SRN model. + * This modifier leads to the {{{CellData}}} cell property being updated at each timestep to deal with Delta-Notch signalling. + */ +#include "DeltaNotchTrackingModifier.hpp" + +/* Having included all the necessary header files, we proceed by defining the test class. + */ +class TestOptogenetics : public AbstractCellBasedTestSuite +{ +public: + + void TestVertexBasedMonolayerWithDeltaNotch() throw (Exception) + { + /* We include the next line because Vertex simulations cannot be run in parallel */ + EXIT_IF_PARALLEL; + + /* First we create a regular vertex mesh. */ + HoneycombVertexMeshGenerator generator(10, 10); + MutableVertexMesh<2,2>* p_mesh = generator.GetMesh(); + + std::vector cells; + MAKE_PTR(WildTypeCellMutationState, p_state); + MAKE_PTR(DifferentiatedCellProliferativeType, p_diff_type); + + for (unsigned elem_index=0; elem_indexGetNumElements(); elem_index++) + { + NoCellCycleModel* p_cc_model = new NoCellCycleModel(); + p_cc_model->SetDimension(2); + + /* We choose to initialise the concentrations to random levels in each cell. */ + std::vector initial_conditions; + initial_conditions.push_back(RandomNumberGenerator::Instance()->ranf()); + initial_conditions.push_back(RandomNumberGenerator::Instance()->ranf()); + MatteoSrnModel* p_srn_model = new MatteoSrnModel(); + p_srn_model->SetInitialConditions(initial_conditions); + + CellPtr p_cell(new Cell(p_state, p_cc_model, p_srn_model)); + p_cell->SetCellProliferativeType(p_diff_type); + double birth_time = -RandomNumberGenerator::Instance()->ranf()*12.0; + p_cell->SetBirthTime(birth_time); + cells.push_back(p_cell); + } + + /* Using the vertex mesh and cells, we create a cell-based population object, and specify which results to + * output to file. */ + VertexBasedCellPopulation<2> cell_population(*p_mesh, cells); + + /* We are now in a position to create and configure the cell-based simulation object, pass a force law to it, + * and run the simulation. We can make the simulation run for longer to see more patterning by increasing the end time. */ + OffLatticeSimulation<2> simulator(cell_population); + simulator.SetOutputDirectory("TestOptogenetics"); + simulator.SetSamplingTimestepMultiple(10); + simulator.SetEndTime(11.0); + + /* Then, we define the modifier class, which automatically updates the values of Delta and Notch within the cells in {{{CellData}}} and passes it to the simulation.*/ + MAKE_PTR(DeltaNotchTrackingModifier<2>, p_modifier); + simulator.AddSimulationModifier(p_modifier); + + MAKE_PTR(MatteoForce<2>, p_force); + simulator.AddForce(p_force); + + /* This modifier assigns target areas to each cell, which are required by the {{{MatteoHondaForce}}}. + */ + MAKE_PTR(ConstantTargetAreaModifier<2>, p_growth_modifier); + simulator.AddSimulationModifier(p_growth_modifier); + simulator.Solve(); + } +}; + +#endif /*TESTOPTOGENETICS_HPP_*/ diff --git a/test/Testmatteo.hpp b/test/Testmatteo.hpp new file mode 100644 index 0000000..ef703c7 --- /dev/null +++ b/test/Testmatteo.hpp @@ -0,0 +1,114 @@ +#ifndef TESTRUNNINGDIFFERENTIALADHESIONSIMULATIONSTUTORIAL_HPP_ +#define TESTRUNNINGDIFFERENTIALADHESIONSIMULATIONSTUTORIAL_HPP_ + +#include +#include "CheckpointArchiveTypes.hpp" +#include "AbstractCellBasedTestSuite.hpp" +#include "HoneycombVertexMeshGenerator.hpp" +#include "ToroidalHoneycombVertexMeshGenerator.hpp" +#include "Toroidal2dVertexMesh.hpp" +#include "CellsGenerator.hpp" +#include "MatteoCellCycleModel.hpp" +#include "CellLabel.hpp" +#include "DifferentiatedCellProliferativeType.hpp" +#include "VertexBasedCellPopulation.hpp" +#include "OffLatticeSimulation.hpp" +#include "SmartPointers.hpp" +#include "FakePetscSetup.hpp" +#include "FarhadifarForce.hpp" +#include "ConstantTargetAreaModifier.hpp" +#include "MatteoModifier.hpp" +#include "CellLabelWriter.hpp" + +class Testmatteo : public AbstractCellBasedTestSuite +{ +public: + + void TestVertexBasedDifferentialAdhesionSimulation() throw (Exception) + { + /* First we create a regular vertex mesh. Here we choose to set the value of the cell rearrangement threshold. */ + HoneycombVertexMeshGenerator generator(20, 20); + MutableVertexMesh<2,2>* p_mesh = generator.GetMesh(); + //ToroidalHoneycombVertexMeshGenerator generator(5, 5); + //Toroidal2dVertexMesh* p_mesh = generator.GetToroidalMesh(); + p_mesh->SetCellRearrangementThreshold(0.1); + + /* We then create some cells using the helper class {{{CellsGenerator}}}. Note that in this simulation + * the cells are all differentiated, and thus no cell division occurs; if we wished, we could modify + * the three lines below in a straightforward manner to incorporate cell proliferation and investigate + * the effect of this on the cell sorting process. */ + std::vector cells; + MAKE_PTR(DifferentiatedCellProliferativeType, p_diff_type); + CellsGenerator cells_generator; + cells_generator.GenerateBasic(cells, p_mesh->GetNumElements(), std::vector(), p_diff_type); + + for(unsigned i=0; iGetCellData()->SetItem("divide", 0); + cells[i]->GetCellData()->SetItem("fitness", 1); + + } + + + /* Using the vertex mesh and cells, we create a cell-based population object, and specify which results to + * output to file. */ + VertexBasedCellPopulation<2> cell_population(*p_mesh, cells); + cell_population.AddCellWriter(); + + /* We randomly label some cells using the cell property {{{CellLabel}}}. We begin by creating a shared pointer to + * this cell property using the helper singleton {{{CellPropertyRegistry}}}. We then loop over the cells and label + * each cell independently with probability 0.5. Note that since the cells have been passed to the + * {{{VertexBasedCellPopulation}}} object, the vector {{{cells}}} above is now empty, so we must use the + * {{{Iterator}}} to loop over cells. */ + boost::shared_ptr p_label(CellPropertyRegistry::Instance()->Get()); + for (AbstractCellPopulation<2>::Iterator cell_iter = cell_population.Begin(); + cell_iter != cell_population.End(); + ++cell_iter) + { + if (RandomNumberGenerator::Instance()->ranf() < 0.1) + { + cell_iter->AddCellProperty(p_label); + } + } + + /* We are now in a position to create and configure the cell-based simulation object. + * We can make the simulation run for longer to see more cell sorting by increasing the end time. */ + OffLatticeSimulation<2> simulator(cell_population); + simulator.SetOutputDirectory("TestMatteo"); + simulator.SetSamplingTimestepMultiple(10); + simulator.SetEndTime(100.0); + + /* Next we create the differential adhesion force law. This builds upon the model of Nagai, Honda and co-workers + * encounted in the TestRunningVertexBasedSimulationsTutorial by allowing different values of the adhesion + * energy parameters depending on the types of two neighbouring cells. Here we interpret the 'type' of a cell + * as whether or not it has the cell property {{{CellLabel}}}; it would be straightforward to create a similar + * force law that took account of a cell's mutation state, for example. Having created the force law, we set the + * values of the parameters. If the adhesion energy for two neighbouring homotypic cells is less than that of two + * heterotypic cells, then we may expect cell sorting to occur, in which the cells of each type will tend to locally + * aggregate over time. */ + MAKE_PTR(FarhadifarForce<2>, p_force); + simulator.AddForce(p_force); + + /* A {{{NagaiHondaForceDifferentialAdhesionForce}}} assumes that each cell has been assigned a target area. + * The {{{ConstantTargetAreaModifier}}} will assign and update the target areas of all cells. + */ + MAKE_PTR(ConstantTargetAreaModifier<2>, p_growth_modifier); + simulator.AddSimulationModifier(p_growth_modifier); + + MAKE_PTR(MatteoModifier<2>, p_modifier); + simulator.AddSimulationModifier(p_modifier); + + /* Finally, we run the simulation. */ + simulator.Solve(); + } + + /* + * EMPTYLINE + * + * To visualize the results, use Paraview. See the UserTutorials/VisualizingWithParaview tutorial for more information. + * + * Load the file {{{/tmp/$USER/testoutput/TestVertexBasedDifferentialAdhesionSimulation/results_from_time_0/results.pvd}}}. + */ +}; + +#endif /*TESTRUNNINGDIFFERENTIALADHESIONSIMULATIONSTUTORIAL_HPP_*/