Chaste  Release::2018.1
LinearSpringWithVariableSpringConstantsForce.cpp
1 /*
2 
3 Copyright (c) 2005-2018, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, are permitted provided that the following conditions are met:
14  * Redistributions of source code must retain the above copyright notice,
15  this list of conditions and the following disclaimer.
16  * Redistributions in binary form must reproduce the above copyright notice,
17  this list of conditions and the following disclaimer in the documentation
18  and/or other materials provided with the distribution.
19  * Neither the name of the University of Oxford nor the names of its
20  contributors may be used to endorse or promote products derived from this
21  software without specific prior written permission.
22 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34 */
35 
36 #include "LinearSpringWithVariableSpringConstantsForce.hpp"
37 #include "MeshBasedCellPopulation.hpp"
38 #include "VanLeeuwen2009WntSwatCellCycleModelHypothesisOne.hpp"
39 #include "VanLeeuwen2009WntSwatCellCycleModelHypothesisTwo.hpp"
40 #include "ApoptoticCellProperty.hpp"
41 #include "BetaCateninOneHitCellMutationState.hpp"
42 #include "ApcTwoHitCellMutationState.hpp"
43 
44 template<unsigned DIM>
47  mUseEdgeBasedSpringConstant(false),
48  mUseMutantSprings(false),
49  mMutantMutantMultiplier(DOUBLE_UNSET),
50  mNormalMutantMultiplier(DOUBLE_UNSET),
51  mUseBCatSprings(false),
52  mUseApoptoticSprings(false),
53  mBetaCatSpringScaler(18.14/6.0), // scale spring constant with beta-catenin level (divided by 6 for heaxagonal cells)
54  mApoptoticSpringTensionStiffness(15.0*0.25),
55  mApoptoticSpringCompressionStiffness(15.0*0.75)
56 {
57 }
58 
59 template<unsigned DIM>
61 {
62 }
63 
64 template<unsigned DIM>
66 {
67  assert(DIM == 2); // LCOV_EXCL_LINE
68  mUseEdgeBasedSpringConstant = useEdgeBasedSpringConstant;
69 }
70 
71 template<unsigned DIM>
72 void LinearSpringWithVariableSpringConstantsForce<DIM>::SetMutantSprings(bool useMutantSprings, double mutantMutantMultiplier, double normalMutantMultiplier)
73 {
74  mUseMutantSprings = useMutantSprings;
75  mMutantMutantMultiplier = mutantMutantMultiplier;
76  mNormalMutantMultiplier = normalMutantMultiplier;
77 }
78 
79 template<unsigned DIM>
81 {
82  mUseBCatSprings = useBCatSprings;
83 }
84 
85 template<unsigned DIM>
87 {
88  mUseApoptoticSprings = useApoptoticSprings;
89 }
90 
91 template<unsigned DIM>
93  unsigned nodeAGlobalIndex,
94  unsigned nodeBGlobalIndex,
95  AbstractCellPopulation<DIM>& rCellPopulation,
96  bool isCloserThanRestLength)
97 {
98 
99  double multiplication_factor = GeneralisedLinearSpringForce<DIM>::VariableSpringConstantMultiplicationFactor(nodeAGlobalIndex,
100  nodeBGlobalIndex,
101  rCellPopulation,
102  isCloserThanRestLength);
103 
104  CellPtr p_cell_A = rCellPopulation.GetCellUsingLocationIndex(nodeAGlobalIndex);
105  CellPtr p_cell_B = rCellPopulation.GetCellUsingLocationIndex(nodeBGlobalIndex);
106 
107  /*
108  * The next code block computes the edge-dependent spring constant as given by equation
109  * (3) in the following reference: van Leeuwen et al. 2009. An integrative computational model
110  * for intestinal tissue renewal. Cell Prolif. 42(5):617-636. doi:10.1111/j.1365-2184.2009.00627.x
111  */
112  if (mUseEdgeBasedSpringConstant)
113  {
114  assert(bool(dynamic_cast<MeshBasedCellPopulation<DIM>*>(&rCellPopulation)));
115  assert(!mUseBCatSprings); // don't want to do both (both account for edge length)
116 
117  multiplication_factor = (static_cast<MeshBasedCellPopulation<DIM>*>(&rCellPopulation))->GetVoronoiEdgeLength(nodeAGlobalIndex, nodeBGlobalIndex)*sqrt(3.0);
118  }
119 
120  if (mUseMutantSprings)
121  {
122  unsigned number_of_mutants = 0;
123 
124  if (p_cell_A->GetMutationState()->IsType<ApcTwoHitCellMutationState>() || p_cell_A->GetMutationState()->IsType<BetaCateninOneHitCellMutationState>())
125  {
126  // If cell A is mutant
127  number_of_mutants++;
128  }
129 
130  if (p_cell_B->GetMutationState()->IsType<ApcTwoHitCellMutationState>() || p_cell_B->GetMutationState()->IsType<BetaCateninOneHitCellMutationState>())
131  {
132  // If cell B is mutant
133  number_of_mutants++;
134  }
135 
136  switch (number_of_mutants)
137  {
138  case 1u:
139  {
140  multiplication_factor *= mNormalMutantMultiplier;
141  break;
142  }
143  case 2u:
144  {
145  multiplication_factor *= mMutantMutantMultiplier;
146  break;
147  }
148  }
149  }
150 
151  /*
152  * The next code block computes the beta-catenin dependent spring constant as given by equation
153  * (4) in the following reference: van Leeuwen et al. 2009. An integrative computational model
154  * for intestinal tissue renewal. Cell Prolif. 42(5):617-636. doi:10.1111/j.1365-2184.2009.00627.x
155  */
156  if (mUseBCatSprings)
157  {
158  assert(bool(dynamic_cast<MeshBasedCellPopulation<DIM>*>(&rCellPopulation)));
159 
160  // If using beta-cat dependent springs, both cell-cycle models had better be VanLeeuwen2009WntSwatCellCycleModel
161  AbstractVanLeeuwen2009WntSwatCellCycleModel* p_model_A = dynamic_cast<AbstractVanLeeuwen2009WntSwatCellCycleModel*>(p_cell_A->GetCellCycleModel());
162  AbstractVanLeeuwen2009WntSwatCellCycleModel* p_model_B = dynamic_cast<AbstractVanLeeuwen2009WntSwatCellCycleModel*>(p_cell_B->GetCellCycleModel());
163 
164  assert(!mUseEdgeBasedSpringConstant); // This already adapts for edge lengths - don't want to do it twice.
165  double beta_cat_cell_1 = p_model_A->GetMembraneBoundBetaCateninLevel();
166  double beta_cat_cell_2 = p_model_B->GetMembraneBoundBetaCateninLevel();
167 
168  MeshBasedCellPopulation<DIM>* p_static_cast_cell_population = (static_cast<MeshBasedCellPopulation<DIM>*>(&rCellPopulation));
169 
170  double perim_cell_1 = p_static_cast_cell_population->GetSurfaceAreaOfVoronoiElement(nodeAGlobalIndex);
171  double perim_cell_2 = p_static_cast_cell_population->GetSurfaceAreaOfVoronoiElement(nodeBGlobalIndex);
172  double edge_length_between_1_and_2 = p_static_cast_cell_population->GetVoronoiEdgeLength(nodeAGlobalIndex, nodeBGlobalIndex);
173 
174  double beta_cat_on_cell_1_edge = beta_cat_cell_1 * edge_length_between_1_and_2 / perim_cell_1;
175  double beta_cat_on_cell_2_edge = beta_cat_cell_2 * edge_length_between_1_and_2 / perim_cell_2;
176 
177  double min_beta_Cat_of_two_cells = std::min(beta_cat_on_cell_1_edge, beta_cat_on_cell_2_edge);
178 
179  multiplication_factor *= min_beta_Cat_of_two_cells / mBetaCatSpringScaler;
180  }
181 
182  if (mUseApoptoticSprings)
183  {
184  bool cell_A_is_apoptotic = p_cell_A->HasCellProperty<ApoptoticCellProperty>();
185  bool cell_B_is_apoptotic = p_cell_B->HasCellProperty<ApoptoticCellProperty>();
186 
187  if (cell_A_is_apoptotic || cell_B_is_apoptotic)
188  {
189  double spring_a_stiffness = 2.0 * this->GetMeinekeSpringStiffness();
190  double spring_b_stiffness = 2.0 * this->GetMeinekeSpringStiffness();
191 
192  if (cell_A_is_apoptotic)
193  {
194  if (!isCloserThanRestLength) // if under tension
195  {
196  spring_a_stiffness = mApoptoticSpringTensionStiffness;
197  }
198  else // if under compression
199  {
200  spring_a_stiffness = mApoptoticSpringCompressionStiffness;
201  }
202  }
203  if (cell_B_is_apoptotic)
204  {
205  if (!isCloserThanRestLength) // if under tension
206  {
207  spring_b_stiffness = mApoptoticSpringTensionStiffness;
208  }
209  else // if under compression
210  {
211  spring_b_stiffness = mApoptoticSpringCompressionStiffness;
212  }
213  }
214 
215  multiplication_factor /= (1.0/spring_a_stiffness + 1.0/spring_b_stiffness)*this->GetMeinekeSpringStiffness();
216  }
217  }
218 
219  return multiplication_factor;
220 }
221 
222 template<unsigned DIM>
224 {
225  // Throw an exception message if not using a MeshBasedCellPopulation
226  if (dynamic_cast<MeshBasedCellPopulation<DIM>*>(&rCellPopulation) == nullptr)
227  {
228  EXCEPTION("LinearSpringWithVariableSpringConstantsForce is to be used with a subclass of MeshBasedCellPopulation only");
229  }
230 
231  MeshBasedCellPopulation<DIM>* p_static_cast_cell_population = static_cast<MeshBasedCellPopulation<DIM>*>(&rCellPopulation);
232 
233  for (typename MeshBasedCellPopulation<DIM>::SpringIterator spring_iterator = p_static_cast_cell_population->SpringsBegin();
234  spring_iterator != p_static_cast_cell_population->SpringsEnd();
235  ++spring_iterator)
236  {
237  unsigned nodeA_global_index = spring_iterator.GetNodeA()->GetIndex();
238  unsigned nodeB_global_index = spring_iterator.GetNodeB()->GetIndex();
239 
240  c_vector<double, DIM> force = this->CalculateForceBetweenNodes(nodeA_global_index, nodeB_global_index, rCellPopulation);
241  c_vector<double, DIM> negative_force = -1.0*force;
242 
243  spring_iterator.GetNodeB()->AddAppliedForceContribution(negative_force);
244  spring_iterator.GetNodeA()->AddAppliedForceContribution(force);
245  }
246 }
247 
248 template<unsigned DIM>
250 {
251  return mBetaCatSpringScaler;
252 }
253 
254 template<unsigned DIM>
256 {
257  assert(betaCatSpringScaler > 0.0);
258  mBetaCatSpringScaler = betaCatSpringScaler;
259 }
260 
261 template<unsigned DIM>
263 {
264  return mApoptoticSpringTensionStiffness;
265 }
266 
267 template<unsigned DIM>
269 {
270  assert(apoptoticSpringTensionStiffness >= 0.0);
271  mApoptoticSpringTensionStiffness = apoptoticSpringTensionStiffness;
272 }
273 
274 template<unsigned DIM>
276 {
277  return mApoptoticSpringCompressionStiffness;
278 }
279 
280 template<unsigned DIM>
282 {
283  assert(apoptoticSpringCompressionStiffness >= 0.0);
284  mApoptoticSpringCompressionStiffness = apoptoticSpringCompressionStiffness;
285 }
286 
287 template<unsigned DIM>
289 {
290  *rParamsFile << "\t\t\t<UseEdgeBasedSpringConstant>" << mUseEdgeBasedSpringConstant << "</UseEdgeBasedSpringConstant>\n";
291  *rParamsFile << "\t\t\t<UseMutantSprings>" << mUseMutantSprings << "</UseMutantSprings>\n";
292  *rParamsFile << "\t\t\t<MutantMutantMultiplier>" << mMutantMutantMultiplier << "</MutantMutantMultiplier>\n";
293  *rParamsFile << "\t\t\t<NormalMutantMultiplier>" << mNormalMutantMultiplier << "</NormalMutantMultiplier>\n";
294  *rParamsFile << "\t\t\t<UseBCatSprings>" << mUseBCatSprings << "</UseBCatSprings>\n";
295  *rParamsFile << "\t\t\t<UseApoptoticSprings>" << mUseApoptoticSprings << "</UseApoptoticSprings>\n";
296  *rParamsFile << "\t\t\t<BetaCatSpringScaler>" << mBetaCatSpringScaler << "</BetaCatSpringScaler>\n";
297  *rParamsFile << "\t\t\t<ApoptoticSpringTensionStiffness>" << mApoptoticSpringTensionStiffness << "</ApoptoticSpringTensionStiffness>\n";
298  *rParamsFile << "\t\t\t<ApoptoticSpringCompressionStiffness>" << mApoptoticSpringCompressionStiffness << "</ApoptoticSpringCompressionStiffness>\n";
299 
300  // Call method on direct parent class
302 }
303 
304 // Explicit instantiation
308 
309 // Serialization for Boost >= 1.36
double GetSurfaceAreaOfVoronoiElement(unsigned index)
double GetVoronoiEdgeLength(unsigned index1, unsigned index2)
void SetApoptoticSpringTensionStiffness(double apoptoticSpringTensionStiffness)
virtual CellPtr GetCellUsingLocationIndex(unsigned index)
double VariableSpringConstantMultiplicationFactor(unsigned nodeAGlobalIndex, unsigned nodeBGlobalIndex, AbstractCellPopulation< DIM > &rCellPopulation, bool isCloserThanRestLength)
#define EXCEPTION(message)
Definition: Exception.hpp:143
const double DOUBLE_UNSET
Definition: Exception.hpp:56
#define EXPORT_TEMPLATE_CLASS_SAME_DIMS(CLASS)
void AddForceContribution(AbstractCellPopulation< DIM > &rCellPopulation)
virtual double VariableSpringConstantMultiplicationFactor(unsigned nodeAGlobalIndex, unsigned nodeBGlobalIndex, AbstractCellPopulation< ELEMENT_DIM, SPACE_DIM > &rCellPopulation, bool isCloserThanRestLength)
void SetMutantSprings(bool useMutantSprings, double mutantMutantMultiplier=2, double normalMutantMultiplier=1.5)
virtual void OutputForceParameters(out_stream &rParamsFile)
void SetApoptoticSpringCompressionStiffness(double apoptoticSpringCompressionStiffness)