CryptProjectionForce.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 
00029 #include "CryptProjectionForce.hpp"
00030 #include "MeshBasedCellPopulation.hpp"
00031 #include "WntConcentration.hpp"
00032 
00033 CryptProjectionForce::CryptProjectionForce()
00034     : GeneralisedLinearSpringForce<2>(),
00035       mIncludeWntChemotaxis(false),
00036       mWntChemotaxisStrength(100.0)
00037 {
00038     mA = WntConcentration<2>::Instance()->GetCryptProjectionParameterA();
00039     mB = WntConcentration<2>::Instance()->GetCryptProjectionParameterB();
00040 }
00041 
00042 CryptProjectionForce::~CryptProjectionForce()
00043 {
00044 }
00045 
00046 void CryptProjectionForce::UpdateNode3dLocationMap(AbstractCellPopulation<2>& rCellPopulation)
00047 {
00048     mNode3dLocationMap.clear();
00049 
00050     c_vector<double, 2> node_location_2d;
00051     c_vector<double, 3> node_location_3d;
00052 
00053     // Only consider nodes corresponding to real cells
00054     for (AbstractCellPopulation<2>::Iterator cell_iter = rCellPopulation.Begin();
00055          cell_iter != rCellPopulation.End();
00056          ++cell_iter)
00057     {
00058         // Get node index
00059         unsigned node_index = rCellPopulation.GetLocationIndexUsingCell(*cell_iter);
00060 
00061         // Get 3D location
00062         node_location_2d = rCellPopulation.GetLocationOfCellCentre(*cell_iter);
00063 
00064         node_location_3d[0] = node_location_2d[0];
00065         node_location_3d[1] = node_location_2d[1];
00066         node_location_3d[2] = CalculateCryptSurfaceHeightAtPoint(node_location_2d);
00067 
00068         // Add to map
00069         mNode3dLocationMap[node_index] = node_location_3d;
00070     }
00071 }
00072 
00073 double CryptProjectionForce::GetA() const
00074 {
00075     return mA;
00076 }
00077 
00078 double CryptProjectionForce::GetB() const
00079 {
00080     return mB;
00081 }
00082 
00083 void CryptProjectionForce::SetWntChemotaxisStrength(double wntChemotaxisStrength)
00084 {
00085     assert(wntChemotaxisStrength >= 0.0);
00086     mWntChemotaxisStrength = wntChemotaxisStrength;
00087 }
00088 double CryptProjectionForce::GetWntChemotaxisStrength()
00089 {
00090     return mWntChemotaxisStrength;
00091 }
00092 
00093 void CryptProjectionForce::SetWntChemotaxis(bool includeWntChemotaxis)
00094 {
00095     mIncludeWntChemotaxis = includeWntChemotaxis;
00096 }
00097 
00098 double CryptProjectionForce::CalculateCryptSurfaceHeightAtPoint(const c_vector<double,2>& rNodeLocation)
00099 {
00100     return mA*pow(norm_2(rNodeLocation), mB); // =z_coord;
00101 }
00102 
00103 double CryptProjectionForce::CalculateCryptSurfaceDerivativeAtPoint(const c_vector<double,2>& rNodeLocation)
00104 {
00105     return mA*mB*pow(norm_2(rNodeLocation), (mB - 1.0));
00106 }
00107 
00108 c_vector<double,2> CryptProjectionForce::CalculateForceBetweenNodes(unsigned nodeAGlobalIndex, unsigned nodeBGlobalIndex, AbstractCellPopulation<2>& rCellPopulation)
00109 {
00110     MeshBasedCellPopulation<2>* p_static_cast_cell_population = static_cast<MeshBasedCellPopulation<2>*>(&rCellPopulation);
00111 
00112     // We should only ever calculate the force between two distinct nodes
00113     assert(nodeAGlobalIndex != nodeBGlobalIndex);
00114 
00115     // Get the node locations in 2D
00116     c_vector<double,2> node_a_location_2d = rCellPopulation.GetNode(nodeAGlobalIndex)->rGetLocation();
00117     c_vector<double,2> node_b_location_2d = rCellPopulation.GetNode(nodeBGlobalIndex)->rGetLocation();
00118 
00119     // "Get the unit vector parallel to the line joining the two nodes" [GeneralisedLinearSpringForce]
00120 
00121     // Create a unit vector in the direction of the 3D spring
00122     c_vector<double,3> unit_difference = mNode3dLocationMap[nodeBGlobalIndex] - mNode3dLocationMap[nodeAGlobalIndex];
00123 
00124     // Calculate the distance between the two nodes
00125     double distance_between_nodes = norm_2(unit_difference);
00126     assert(distance_between_nodes > 0);
00127     assert(!std::isnan(distance_between_nodes));
00128 
00129     unit_difference /= distance_between_nodes;
00130 
00131     // If mUseCutOffLength has been set, then there is zero force between
00132     // two nodes located a distance apart greater than mUseCutOffLength
00133     if (this->mUseCutOffLength)
00134     {
00135         if (distance_between_nodes >= mMechanicsCutOffLength)
00136         {
00137             // Return zero (2D projected) force
00138             return zero_vector<double>(2);
00139         }
00140     }
00141 
00142     // Calculate the rest length of the spring connecting the two nodes
00143 
00144     double rest_length = 1.0;
00145 
00146     CellPtr p_cell_A = rCellPopulation.GetCellUsingLocationIndex(nodeAGlobalIndex);
00147     CellPtr p_cell_B = rCellPopulation.GetCellUsingLocationIndex(nodeBGlobalIndex);
00148 
00149     double ageA = p_cell_A->GetAge();
00150     double ageB = p_cell_B->GetAge();
00151 
00152     assert(!std::isnan(ageA));
00153     assert(!std::isnan(ageB));
00154 
00155     /*
00156      * If the cells are both newly divided, then the rest length of the spring
00157      * connecting them grows linearly with time, until mMeinekeSpringGrowthDuration hour after division.
00158      */
00159     if (ageA < mMeinekeSpringGrowthDuration && ageB < mMeinekeSpringGrowthDuration)
00160     {
00161         /*
00162          * The spring rest length increases from a predefined small parameter
00163          * to a normal rest length of 1.0, over a period of one hour.
00164          */
00165         std::pair<CellPtr,CellPtr> cell_pair = p_static_cast_cell_population->CreateCellPair(p_cell_A, p_cell_B);
00166         if (p_static_cast_cell_population->IsMarkedSpring(cell_pair))
00167         {
00168             double lambda = mMeinekeDivisionRestingSpringLength;
00169             rest_length = lambda + (1.0 - lambda) * ageA/mMeinekeSpringGrowthDuration;
00170         }
00171         if (ageA+SimulationTime::Instance()->GetTimeStep() >= mMeinekeSpringGrowthDuration)
00172         {
00173             // This spring is about to go out of scope
00174             p_static_cast_cell_population->UnmarkSpring(cell_pair);
00175         }
00176     }
00177 
00178     double a_rest_length = rest_length*0.5;
00179     double b_rest_length = a_rest_length;
00180 
00181     /*
00182      * If either of the cells has begun apoptosis, then the length of the spring
00183      * connecting them decreases linearly with time.
00184      */
00185     if (p_cell_A->HasApoptosisBegun())
00186     {
00187         double time_until_death_a = p_cell_A->GetTimeUntilDeath();
00188         a_rest_length = a_rest_length * time_until_death_a / p_cell_A->GetApoptosisTime();
00189     }
00190     if (p_cell_B->HasApoptosisBegun())
00191     {
00192         double time_until_death_b = p_cell_B->GetTimeUntilDeath();
00193         b_rest_length = b_rest_length * time_until_death_b / p_cell_B->GetApoptosisTime();
00194     }
00195 
00196     rest_length = a_rest_length + b_rest_length;
00197 
00198     // Assert that the rest length does not exceed 1
00199     assert(rest_length <= 1.0+1e-12);
00200 
00201     bool is_closer_than_rest_length = true;
00202 
00203     if (distance_between_nodes - rest_length > 0)
00204     {
00205         is_closer_than_rest_length = false;
00206     }
00207 
00208     /*
00209      * Although in this class the 'spring constant' is a constant parameter, in
00210      * subclasses it can depend on properties of each of the cells.
00211      */
00212     double multiplication_factor = 1.0;
00213     multiplication_factor *= VariableSpringConstantMultiplicationFactor(nodeAGlobalIndex, nodeBGlobalIndex, rCellPopulation, is_closer_than_rest_length);
00214 
00215     // Calculate the 3D force between the two points
00216     c_vector<double,3> force_between_nodes = multiplication_factor * this->GetMeinekeSpringStiffness() * unit_difference * (distance_between_nodes - rest_length);
00217 
00218     // Calculate an outward normal unit vector to the tangent plane of the crypt surface at the 3D point corresponding to node B
00219     c_vector<double,3> outward_normal_unit_vector;
00220 
00221     double dfdr = CalculateCryptSurfaceDerivativeAtPoint(node_b_location_2d);
00222     double theta_B = atan2(node_b_location_2d[1], node_b_location_2d[0]); // use atan2 to determine the quadrant
00223     double normalization_factor = sqrt(1 + dfdr*dfdr);
00224 
00225     outward_normal_unit_vector[0] = dfdr*cos(theta_B)/normalization_factor;
00226     outward_normal_unit_vector[1] = dfdr*sin(theta_B)/normalization_factor;
00227     outward_normal_unit_vector[2] = -1.0/normalization_factor;
00228 
00229     // Calculate the projection of the force onto the plane z=0
00230     c_vector<double,2> projected_force_between_nodes_2d;
00231     double force_dot_normal = inner_prod(force_between_nodes, outward_normal_unit_vector);
00232 
00233     for (unsigned i=0; i<2; i++)
00234     {
00235         projected_force_between_nodes_2d[i] = force_between_nodes[i]
00236                                               - force_dot_normal*outward_normal_unit_vector[i]
00237                                               + force_dot_normal*outward_normal_unit_vector[2];
00238     }
00239 
00240     return projected_force_between_nodes_2d;
00241 }
00242 
00243 void CryptProjectionForce::AddForceContribution(std::vector<c_vector<double,2> >& rForces,
00244                                                 AbstractCellPopulation<2>& rCellPopulation)
00245 {
00246     // First work out the 3D location of each cell
00247     UpdateNode3dLocationMap(rCellPopulation);
00248 
00249     // Throw an exception message if not using a MeshBasedCellPopulation
00250     if (dynamic_cast<MeshBasedCellPopulation<2>*>(&rCellPopulation) == NULL)
00251     {
00252         EXCEPTION("CryptProjectionForce is to be used with a subclass of MeshBasedCellPopulation only");
00253     }
00254 
00255     MeshBasedCellPopulation<2>* p_static_cast_cell_population = static_cast<MeshBasedCellPopulation<2>*>(&rCellPopulation);
00256 
00257     for (MeshBasedCellPopulation<2>::SpringIterator spring_iterator = p_static_cast_cell_population->SpringsBegin();
00258          spring_iterator != p_static_cast_cell_population->SpringsEnd();
00259          ++spring_iterator)
00260     {
00261         unsigned nodeA_global_index = spring_iterator.GetNodeA()->GetIndex();
00262         unsigned nodeB_global_index = spring_iterator.GetNodeB()->GetIndex();
00263 
00264         c_vector<double, 2> force = CalculateForceBetweenNodes(nodeA_global_index, nodeB_global_index, rCellPopulation);
00265 
00266         rForces[nodeB_global_index] -= force;
00267         rForces[nodeA_global_index] += force;
00268     }
00269 
00270     if (mIncludeWntChemotaxis)
00271     {
00272         assert(WntConcentration<2>::Instance()->IsWntSetUp());
00273 
00274         for (AbstractCellPopulation<2>::Iterator cell_iter = rCellPopulation.Begin();
00275              cell_iter != rCellPopulation.End();
00276              ++cell_iter)
00277         {
00278             if (cell_iter->GetCellCycleModel()->GetCellProliferativeType()==STEM)
00279             {
00280                 c_vector<double, 2> wnt_chemotactic_force = mWntChemotaxisStrength*WntConcentration<2>::Instance()->GetWntGradient(*cell_iter);
00281                 unsigned index = rCellPopulation.GetLocationIndexUsingCell(*cell_iter);
00282 
00283                 rForces[index] += wnt_chemotactic_force;
00284             }
00285         }
00286     }
00287 }
00288 
00289 void CryptProjectionForce::OutputForceParameters(out_stream& rParamsFile)
00290 {
00291     *rParamsFile << "\t\t\t<A>" << mA << "</A> \n";
00292     *rParamsFile << "\t\t\t<B>" << mB << "</B> \n";
00293     *rParamsFile << "\t\t\t<IncludeWntChemotaxis>" << mIncludeWntChemotaxis << "</IncludeWntChemotaxis> \n";
00294     *rParamsFile << "\t\t\t<WntChemotaxisStrength>" << mWntChemotaxisStrength << "</WntChemotaxisStrength> \n";
00295 
00296     // Call method on direct parent class
00297     GeneralisedLinearSpringForce<2>::OutputForceParameters(rParamsFile);
00298 }
00299 
00300 // Serialization for Boost >= 1.36
00301 #include "SerializationExportWrapperForCpp.hpp"
00302 CHASTE_CLASS_EXPORT(CryptProjectionForce)

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