NagaiHondaForce.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 #include "NagaiHondaForce.hpp"
00029 
00030 template<unsigned DIM>
00031 NagaiHondaForce<DIM>::NagaiHondaForce()
00032    : AbstractForce<DIM>(),
00033      mNagaiHondaDeformationEnergyParameter(100.0), // This is 1.0 in the Nagai & Honda paper
00034      mNagaiHondaMembraneSurfaceEnergyParameter(10.0), // This is 0.1 the Nagai & Honda paper
00035      mNagaiHondaCellCellAdhesionEnergyParameter(1.0), // This is 0.01 the Nagai & Honda paper
00036      mNagaiHondaCellBoundaryAdhesionEnergyParameter(1.0), // This is 0.01 the Nagai & Honda paper
00037      mMatureCellTargetArea(1.0)
00038 {
00039 }
00040 
00041 template<unsigned DIM>
00042 NagaiHondaForce<DIM>::~NagaiHondaForce()
00043 {
00044 }
00045 
00046 template<unsigned DIM>
00047 void NagaiHondaForce<DIM>::AddForceContribution(std::vector<c_vector<double, DIM> >& rForces,
00048                                                 AbstractCellPopulation<DIM>& rCellPopulation)
00049 {
00050     // Throw an exception message if not using a VertexBasedCellPopulation
00051     if (dynamic_cast<VertexBasedCellPopulation<DIM>*>(&rCellPopulation) == NULL)
00052     {
00053         EXCEPTION("NagaiHondaForce is to be used with a VertexBasedCellPopulation only");
00054     }
00055 
00056     // Helper variable that is a static cast of the cell population
00057     VertexBasedCellPopulation<DIM>* p_cell_population = static_cast<VertexBasedCellPopulation<DIM>*>(&rCellPopulation);
00058 
00059     // Iterate over vertices in the cell population
00060     for (unsigned node_index=0; node_index<p_cell_population->GetNumNodes(); node_index++)
00061     {
00062         // Compute the force on this node
00063 
00064         /*
00065          * The force on this Node is given by the gradient of the total free
00066          * energy of the CellPopulation, evaluated at the position of the vertex. This
00067          * free energy is the sum of the free energies of all CellPtrs in
00068          * the cell population. The free energy of each CellPtr is comprised of three
00069          * parts - a cell deformation energy, a membrane surface tension energy
00070          * and an adhesion energy.
00071          *
00072          * Note that since the movement of this Node only affects the free energy
00073          * of the three CellPtrs containing it, we can just consider the
00074          * contributions to the free energy gradient from each of these three
00075          * CellPtrs.
00076          */
00077 
00078         c_vector<double, DIM> deformation_contribution = zero_vector<double>(DIM);
00079         c_vector<double, DIM> membrane_surface_tension_contribution = zero_vector<double>(DIM);
00080         c_vector<double, DIM> adhesion_contribution = zero_vector<double>(DIM);
00081 
00082         // Find the indices of the elements owned by this node
00083         std::set<unsigned> containing_elem_indices = p_cell_population->GetNode(node_index)->rGetContainingElementIndices();
00084 
00085         // Iterate over these elements
00086         for (std::set<unsigned>::iterator iter = containing_elem_indices.begin();
00087              iter != containing_elem_indices.end();
00088              ++iter)
00089         {
00090             // Get this element and its index
00091             VertexElement<DIM, DIM>* p_element = p_cell_population->GetElement(*iter);
00092             unsigned element_index = p_element->GetIndex();
00093 
00094             // Find the local index of this node in this element
00095             unsigned local_index = p_element->GetNodeLocalIndex(node_index);
00096 
00097             /******** Start of deformation force calculation ********/
00098 
00099             // Compute the area of this element and its gradient at this node
00100             double element_area = p_cell_population->rGetMesh().GetVolumeOfElement(*iter);
00101             c_vector<double, DIM> element_area_gradient = p_cell_population->rGetMesh().GetAreaGradientOfElementAtNode(p_element, local_index);
00102 
00103             // Get the target area of the cell
00104             double cell_target_area = GetTargetAreaOfCell(p_cell_population->GetCellUsingLocationIndex(element_index));
00105 
00106             // Add the force contribution from this cell's deformation energy (note the minus sign)
00107             deformation_contribution -= 2*GetNagaiHondaDeformationEnergyParameter()*(element_area - cell_target_area)*element_area_gradient;
00108 
00109             /******** End of deformation force calculation *************/
00110 
00111             /******** Start of membrane force calculation ***********/
00112 
00113             // Compute the perimeter of the element and its gradient at this node
00114             double element_perimeter = p_cell_population->rGetMesh().GetSurfaceAreaOfElement(*iter);
00115             c_vector<double, DIM> element_perimeter_gradient = p_cell_population->rGetMesh().GetPerimeterGradientOfElementAtNode(p_element, local_index);
00116 
00117             // Get the target perimeter of the cell
00118             double cell_target_perimeter = 2*sqrt(M_PI*cell_target_area);
00119 
00120             // Add the force contribution from this cell's membrane surface tension (note the minus sign)
00121             membrane_surface_tension_contribution -= 2*GetNagaiHondaMembraneSurfaceEnergyParameter()*(element_perimeter - cell_target_perimeter)*element_perimeter_gradient;
00122 
00123             /******** End of membrane force calculation **********/
00124 
00125             /******** Start of adhesion force calculation ***********/
00126 
00127             // Get the current, previous and next nodes in this element
00128             Node<DIM>* p_current_node = p_element->GetNode(local_index);
00129 
00130             unsigned previous_node_local_index = (p_element->GetNumNodes()+local_index-1)%(p_element->GetNumNodes());
00131             Node<DIM>* p_previous_node = p_element->GetNode(previous_node_local_index);
00132 
00133             unsigned next_node_local_index = (local_index+1)%(p_element->GetNumNodes());
00134             Node<DIM>* p_next_node = p_element->GetNode(next_node_local_index);
00135 
00136             // Compute the adhesion parameter for each of these edges
00137             double previous_edge_adhesion_parameter;
00138             double next_edge_adhesion_parameter;
00139 
00140             previous_edge_adhesion_parameter = GetAdhesionParameter(p_previous_node, p_current_node, *p_cell_population);
00141             next_edge_adhesion_parameter = GetAdhesionParameter(p_current_node, p_next_node, *p_cell_population);
00142 
00143             // Compute the gradient of the edge of the cell ending in this node
00144             c_vector<double, DIM> previous_edge_gradient = p_cell_population->rGetMesh().GetPreviousEdgeGradientOfElementAtNode(p_element, local_index);
00145 
00146             // Compute the gradient of the edge of the cell starting in this node
00147             c_vector<double, DIM> next_edge_gradient = p_cell_population->rGetMesh().GetNextEdgeGradientOfElementAtNode(p_element, local_index);
00148 
00149             // Add the force contribution from cell-cell and cell-boundary adhesion (note the minus sign)
00150             adhesion_contribution -= previous_edge_adhesion_parameter*previous_edge_gradient + next_edge_adhesion_parameter*next_edge_gradient;
00151 
00152             /******** End of adhesion force calculation *************/
00153         }
00154 
00155         c_vector<double, DIM> force_on_node = deformation_contribution +
00156                                               membrane_surface_tension_contribution +
00157                                               adhesion_contribution;
00158 
00159         rForces[node_index] += force_on_node;
00160     }
00161 }
00162 
00163 template<unsigned DIM>
00164 double NagaiHondaForce<DIM>::GetAdhesionParameter(Node<DIM>* pNodeA, Node<DIM>* pNodeB, VertexBasedCellPopulation<DIM>& rVertexCellPopulation)
00165 {
00166     // Find the indices of the elements owned by each node
00167     std::set<unsigned> elements_containing_nodeA = pNodeA->rGetContainingElementIndices();
00168     std::set<unsigned> elements_containing_nodeB = pNodeB->rGetContainingElementIndices();
00169 
00170     // Find common elements
00171     std::set<unsigned> shared_elements;
00172     std::set_intersection(elements_containing_nodeA.begin(),
00173                           elements_containing_nodeA.end(),
00174                           elements_containing_nodeB.begin(),
00175                           elements_containing_nodeB.end(),
00176                           std::inserter(shared_elements, shared_elements.begin()));
00177 
00178     // Check that the nodes have a common edge
00179     assert(!shared_elements.empty());
00180 
00181     double adhesion_parameter = GetNagaiHondaCellCellAdhesionEnergyParameter();
00182 
00183     // If the edge corresponds to a single element, then the cell is on the boundary
00184     if (shared_elements.size() == 1)
00185     {
00186         adhesion_parameter = GetNagaiHondaCellBoundaryAdhesionEnergyParameter();
00187     }
00188 
00189     return adhesion_parameter;
00190 }
00191 
00192 template<unsigned DIM>
00193 double NagaiHondaForce<DIM>::GetNagaiHondaDeformationEnergyParameter()
00194 {
00195     return mNagaiHondaDeformationEnergyParameter;
00196 }
00197 
00198 template<unsigned DIM>
00199 double NagaiHondaForce<DIM>::GetNagaiHondaMembraneSurfaceEnergyParameter()
00200 {
00201     return mNagaiHondaMembraneSurfaceEnergyParameter;
00202 }
00203 
00204 template<unsigned DIM>
00205 double NagaiHondaForce<DIM>::GetNagaiHondaCellCellAdhesionEnergyParameter()
00206 {
00207     return mNagaiHondaCellCellAdhesionEnergyParameter;
00208 }
00209 
00210 template<unsigned DIM>
00211 double NagaiHondaForce<DIM>::GetNagaiHondaCellBoundaryAdhesionEnergyParameter()
00212 {
00213     return mNagaiHondaCellBoundaryAdhesionEnergyParameter;
00214 }
00215 
00216 template<unsigned DIM>
00217 void NagaiHondaForce<DIM>::SetNagaiHondaDeformationEnergyParameter(double deformationEnergyParameter)
00218 {
00219     mNagaiHondaDeformationEnergyParameter = deformationEnergyParameter;
00220 }
00221 
00222 template<unsigned DIM>
00223 void NagaiHondaForce<DIM>::SetNagaiHondaMembraneSurfaceEnergyParameter(double membraneSurfaceEnergyParameter)
00224 {
00225     mNagaiHondaMembraneSurfaceEnergyParameter = membraneSurfaceEnergyParameter;
00226 }
00227 
00228 template<unsigned DIM>
00229 void NagaiHondaForce<DIM>::SetNagaiHondaCellCellAdhesionEnergyParameter(double cellCellAdhesionEnergyParameter)
00230 {
00231     mNagaiHondaCellCellAdhesionEnergyParameter = cellCellAdhesionEnergyParameter;
00232 }
00233 
00234 template<unsigned DIM>
00235 void NagaiHondaForce<DIM>::SetNagaiHondaCellBoundaryAdhesionEnergyParameter(double cellBoundaryAdhesionEnergyParameter)
00236 {
00237     mNagaiHondaCellBoundaryAdhesionEnergyParameter = cellBoundaryAdhesionEnergyParameter;
00238 }
00239 
00240 template<unsigned DIM>
00241 double NagaiHondaForce<DIM>::GetTargetAreaOfCell(const CellPtr pCell) const
00242 {
00243     // Get target area A of a healthy cell in S, G2 or M phase
00244     double cell_target_area = mMatureCellTargetArea;
00245 
00246     double cell_age = pCell->GetAge();
00247     double g1_duration = pCell->GetCellCycleModel()->GetG1Duration();
00248 
00249     // If the cell is differentiated then its G1 duration is infinite
00250     if (g1_duration == DBL_MAX) // don't use magic number, compare to DBL_MAX
00251     {
00252         // This is just for fixed cell-cycle models, need to work out how to find the g1 duration
00253         g1_duration = pCell->GetCellCycleModel()->GetTransitCellG1Duration();
00254     }
00255 
00256     if (pCell->HasCellProperty<ApoptoticCellProperty>())
00257     {
00258         // Age of cell when apoptosis begins
00259         if (pCell->GetStartOfApoptosisTime() - pCell->GetBirthTime() < g1_duration)
00260         {
00261             cell_target_area *= 0.5*(1 + (pCell->GetStartOfApoptosisTime() - pCell->GetBirthTime())/g1_duration);
00262         }
00263 
00264         // The target area of an apoptotic cell decreases linearly to zero (and past it negative)
00265         cell_target_area = cell_target_area - 0.5*cell_target_area/(pCell->GetApoptosisTime())*(SimulationTime::Instance()->GetTime()-pCell->GetStartOfApoptosisTime());
00266 
00267         // Don't allow a negative target area
00268         if (cell_target_area < 0)
00269         {
00270             cell_target_area = 0;
00271         }
00272     }
00273     else
00274     {
00275         // The target area of a proliferating cell increases linearly from A/2 to A over the course of the G1 phase
00276         if (cell_age < g1_duration)
00277         {
00278             cell_target_area *= 0.5*(1 + cell_age/g1_duration);
00279         }
00280     }
00281 
00282     return cell_target_area;
00283 }
00284 
00285 template<unsigned DIM>
00286 double NagaiHondaForce<DIM>::GetMatureCellTargetArea() const
00287 {
00288     return mMatureCellTargetArea;
00289 }
00290 
00291 template<unsigned DIM>
00292 void NagaiHondaForce<DIM>::SetMatureCellTargetArea(double matureCellTargetArea)
00293 {
00294     assert(matureCellTargetArea >= 0.0);
00295     mMatureCellTargetArea = matureCellTargetArea;
00296 }
00297 
00298 template<unsigned DIM>
00299 void NagaiHondaForce<DIM>::OutputForceParameters(out_stream& rParamsFile)
00300 {
00301     *rParamsFile << "\t\t\t<NagaiHondaDeformationEnergyParameter>" << mNagaiHondaDeformationEnergyParameter << "</NagaiHondaDeformationEnergyParameter> \n";
00302     *rParamsFile << "\t\t\t<NagaiHondaMembraneSurfaceEnergyParameter>" << mNagaiHondaMembraneSurfaceEnergyParameter << "</NagaiHondaMembraneSurfaceEnergyParameter> \n";
00303     *rParamsFile << "\t\t\t<NagaiHondaCellCellAdhesionEnergyParameter>" << mNagaiHondaCellCellAdhesionEnergyParameter << "</NagaiHondaCellCellAdhesionEnergyParameter> \n";
00304     *rParamsFile << "\t\t\t<NagaiHondaCellBoundaryAdhesionEnergyParameter>" << mNagaiHondaCellBoundaryAdhesionEnergyParameter << "</NagaiHondaCellBoundaryAdhesionEnergyParameter> \n";
00305     *rParamsFile << "\t\t\t<MatureCellTargetArea>" << mMatureCellTargetArea << "</MatureCellTargetArea> \n";
00306 
00307     // Call method on direct parent class
00308     AbstractForce<DIM>::OutputForceParameters(rParamsFile);
00309 }
00310 
00311 
00313 // Explicit instantiation
00315 
00316 template class NagaiHondaForce<1>;
00317 template class NagaiHondaForce<2>;
00318 template class NagaiHondaForce<3>;
00319 
00320 
00321 // Serialization for Boost >= 1.36
00322 #include "SerializationExportWrapperForCpp.hpp"
00323 EXPORT_TEMPLATE_CLASS_SAME_DIMS(NagaiHondaForce)

Generated on Tue May 31 14:31:40 2011 for Chaste by  doxygen 1.5.5