Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
LinearSpringWithVariableSpringConstantsForce.cpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, 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
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF 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
44template<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
59template<unsigned DIM>
63
64template<unsigned DIM>
66{
67 assert(DIM == 2); // LCOV_EXCL_LINE
68 mUseEdgeBasedSpringConstant = useEdgeBasedSpringConstant;
69}
70
71template<unsigned DIM>
72void LinearSpringWithVariableSpringConstantsForce<DIM>::SetMutantSprings(bool useMutantSprings, double mutantMutantMultiplier, double normalMutantMultiplier)
73{
74 mUseMutantSprings = useMutantSprings;
75 mMutantMutantMultiplier = mutantMutantMultiplier;
76 mNormalMutantMultiplier = normalMutantMultiplier;
77}
78
79template<unsigned DIM>
81{
82 mUseBCatSprings = useBCatSprings;
83}
84
85template<unsigned DIM>
87{
88 mUseApoptoticSprings = useApoptoticSprings;
89}
90
91template<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
222template<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
248template<unsigned DIM>
250{
251 return mBetaCatSpringScaler;
252}
253
254template<unsigned DIM>
256{
257 assert(betaCatSpringScaler > 0.0);
258 mBetaCatSpringScaler = betaCatSpringScaler;
259}
260
261template<unsigned DIM>
263{
264 return mApoptoticSpringTensionStiffness;
265}
266
267template<unsigned DIM>
269{
270 assert(apoptoticSpringTensionStiffness >= 0.0);
271 mApoptoticSpringTensionStiffness = apoptoticSpringTensionStiffness;
272}
273
274template<unsigned DIM>
276{
277 return mApoptoticSpringCompressionStiffness;
278}
279
280template<unsigned DIM>
282{
283 assert(apoptoticSpringCompressionStiffness >= 0.0);
284 mApoptoticSpringCompressionStiffness = apoptoticSpringCompressionStiffness;
285}
286
287template<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
const double DOUBLE_UNSET
Definition Exception.hpp:57
#define EXCEPTION(message)
#define EXPORT_TEMPLATE_CLASS_SAME_DIMS(CLASS)
virtual CellPtr GetCellUsingLocationIndex(unsigned index)
virtual double VariableSpringConstantMultiplicationFactor(unsigned nodeAGlobalIndex, unsigned nodeBGlobalIndex, AbstractCellPopulation< ELEMENT_DIM, SPACE_DIM > &rCellPopulation, bool isCloserThanRestLength)
virtual void OutputForceParameters(out_stream &rParamsFile)
void AddForceContribution(AbstractCellPopulation< DIM > &rCellPopulation)
void SetMutantSprings(bool useMutantSprings, double mutantMutantMultiplier=2, double normalMutantMultiplier=1.5)
double VariableSpringConstantMultiplicationFactor(unsigned nodeAGlobalIndex, unsigned nodeBGlobalIndex, AbstractCellPopulation< DIM > &rCellPopulation, bool isCloserThanRestLength)
void SetApoptoticSpringCompressionStiffness(double apoptoticSpringCompressionStiffness)
void SetApoptoticSpringTensionStiffness(double apoptoticSpringTensionStiffness)
double GetVoronoiEdgeLength(unsigned index1, unsigned index2)
double GetSurfaceAreaOfVoronoiElement(unsigned index)