Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
WntCellCycleOdeSystem.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 "WntCellCycleOdeSystem.hpp"
37#include "CellwiseOdeSystemInformation.hpp"
38
39// These #includes are needed for the constructor and EvaluateYDerivatives()
40#include "ApcOneHitCellMutationState.hpp"
41#include "ApcTwoHitCellMutationState.hpp"
42#include "BetaCateninOneHitCellMutationState.hpp"
43
45 boost::shared_ptr<AbstractCellMutationState> pMutationState,
46 std::vector<double> stateVariables)
48 mpMutationState(pMutationState),
49 mWntLevel(wntLevel)
50{
52
67 Init(); // set up parameter values
68
69 // Set up a Wnt signalling pathway in a steady state
70 double destruction_level = ma5d/(ma4d*wntLevel+ma5d);
71 double beta_cat_level_1 = -1.0;
72 double beta_cat_level_2 = -1.0;
73
74 if (!mpMutationState)
75 {
76 // No mutations specified
77 }
79 {
80 // APC +/- : only half are active
81 beta_cat_level_1 = 0.5*ma2d/(ma2d+0.5*ma3d*destruction_level);
82 beta_cat_level_2 = 0.5*ma2d/(ma2d+0.5*ma3d*destruction_level);
83 }
85 {
86 // APC -/-
87 destruction_level = 0.0; // no active destruction complex
88 beta_cat_level_1 = 0.5; // fully active beta-catenin
89 beta_cat_level_2 = 0.5; // fully active beta-catenin
90 }
92 {
93 // Beta-cat delta 45
94 beta_cat_level_1 = 0.5*ma2d/(ma2d+ma3d*destruction_level);
95 beta_cat_level_2 = 0.5;
96 }
97 else
98 {
99 // healthy cells
100 beta_cat_level_1 = 0.5*ma2d/(ma2d+ma3d*destruction_level);
101 beta_cat_level_2 = 0.5*ma2d/(ma2d+ma3d*destruction_level);
102 }
103
104 // Cell-specific initial conditions
105 SetDefaultInitialCondition(5, destruction_level);
106 SetDefaultInitialCondition(6, beta_cat_level_1);
107 SetDefaultInitialCondition(7, beta_cat_level_2);
108 SetDefaultInitialCondition(8, wntLevel);
109
110 if (stateVariables != std::vector<double>())
111 {
112 SetStateVariables(stateVariables);
113 }
114}
115
116void WntCellCycleOdeSystem::SetMutationState(boost::shared_ptr<AbstractCellMutationState> pMutationState)
117{
118 mpMutationState = pMutationState;
119}
120
125
127{
128 // Initialise model parameter values
129 // Swat (2004) Parameters
130 double k1 = 1.0;
131 double k2 = 1.6;
132 double k3 = 0.05;
133 double k16 = 0.4;
134 double k34 = 0.04;
135 double k43 = 0.01;
136 double k61 = 0.3;
137 double k23 = 0.3;
138 double a = 0.04;
139 double J11 = 0.5;
140 double J12 = 5.0;
141 double J61 = 5.0;
142 double J62 = 8.0;
143 double J13 = 0.002;
144 double J63 = 2.0;
145 double Km1 = 0.5;
146 double Km2 = 4.0;
147 double Km4 = 0.3;
148 double kp = 0.05;
149 double phi_pRb = 0.005;
150 double phi_E2F1 = 0.1;
151 double phi_CycDi = 0.023;
152 double phi_CycDa = 0.03;
153 double phi_pRbp = 0.06;
154
155 // Mirams et al. parameter values
156 double a1 = 0.423;
157 double a2 = 2.57e-4;
158 double a3 = 1.72;
159 double a4 = 10.0;
160 double a5 = 0.5;
161 double WntMax = 10.0;
162 double mitogenic_factorF = 6.0e-4;
163 double APC_Total = 0.02;
164
165 // Non-dimensionalise...
166 mk2d = k2/(Km2*phi_E2F1);
167 mk3d = k3*a1*mitogenic_factorF/(Km4*phi_E2F1*a2);
168 mk34d = k34/phi_E2F1;
169 mk43d = k43/phi_E2F1;
170 mk23d = k23*Km2/(Km4*phi_E2F1);
171 mad = a/Km2;
172 mJ11d = J11*phi_E2F1/k1;
173 mJ12d = J12*phi_E2F1/k1;
174 mJ13d = J13*phi_E2F1/k1;
175 mJ61d = J61*phi_E2F1/k1;
176 mJ62d = J62*phi_E2F1/k1;
177 mJ63d = J63*phi_E2F1/k1;
178 mKm1d = Km1/Km2;
179 mkpd = kp/(Km2*phi_E2F1);
180 mphi_r = phi_pRb/phi_E2F1;
181 mphi_i = phi_CycDi/phi_E2F1;
182 mphi_j = phi_CycDa/phi_E2F1;
183 mphi_p = phi_pRbp/phi_E2F1;
184 ma2d = a2/phi_E2F1;
185 ma3d = a3*APC_Total/phi_E2F1;
186 ma4d = a4*WntMax/phi_E2F1;
187 ma5d = a5/phi_E2F1;
188 mk16d = k16*Km4/phi_E2F1;
189 mk61d = k61/phi_E2F1;
190 mPhiE2F1 = phi_E2F1;
191}
192
193void WntCellCycleOdeSystem::EvaluateYDerivatives(double time, const std::vector<double>& rY, std::vector<double>& rDY)
194{
195 double r = rY[0];
196 double e = rY[1];
197 double i = rY[2];
198 double j = rY[3];
199 double p = rY[4];
200 double c = rY[5];
201 double b1 = rY[6];
202 double b2 = rY[7];
203 double wnt_level = rY[8];
204
205 double dx1 = 0.0;
206 double dx2 = 0.0;
207 double dx3 = 0.0;
208 double dx4 = 0.0;
209 double dx5 = 0.0;
210 double dx6 = 0.0;
211 double dx7 = 0.0;
212 double dx8 = 0.0;
213
214 /*
215 * The variables are
216 * 1. r = pRb
217 * 2. e = E2F1
218 * 3. i = CycD (inactive)
219 * 4. j = CycD (active)
220 * 5. p = pRb-p
221 * 6. c = APC (Active)
222 * 7. b = Beta-Catenin
223 */
224
225 // Bit back-to-front, but work out the Wnt section first...
226
227 // Mutations take effect by altering the level of beta-catenin
228 if (!mpMutationState)
229 {
230 // No mutations specified
231 }
232 else if (mpMutationState->IsType<ApcOneHitCellMutationState>()) // APC +/-
233 {
234 dx6 = ma5d*(1.0-c) - ma4d*wnt_level*c;
235 dx7 = ma2d*(0.5-b1) - 0.5*ma3d*b1*c;
236 dx8 = ma2d*(0.5-b2) - 0.5*ma3d*b2*c;
237 }
238 else if (mpMutationState->IsType<ApcTwoHitCellMutationState>()) // APC -/-
239 {
240 dx6 = 0.0;
241 dx7 = ma2d*(0.5-b1);
242 dx8 = ma2d*(0.5-b2);
243 }
244 else if (mpMutationState->IsType<BetaCateninOneHitCellMutationState>()) // Beta-Cat D45
245 {
246 dx6 = ma5d*(1.0-c) - ma4d*wnt_level*c;
247 dx7 = ma2d*(0.5-b1) - ma3d*b1*c;
248 dx8 = ma2d*(0.5-b2);
249 }
250 else
251 {
252 // da
253 dx6 = ma5d*(1.0-c) - ma4d*wnt_level*c;
254 // db
255 dx7 = ma2d*(0.5-b1) - ma3d*b1*c;
256 dx8 = ma2d*(0.5-b2) - ma3d*b2*c;
257 }
258
259 // Now the cell cycle stuff...
260
261 // dr
262 dx1 = e/(mKm1d+e)*mJ11d/(mJ11d+r)*mJ61d/(mJ61d+p) - mk16d*r*j+mk61d*p-mphi_r*r;
263 // de
264 dx2 = mkpd+mk2d*(mad*mad+e*e)/(1+e*e)*mJ12d/(mJ12d+r)*mJ62d/(mJ62d+p) - e;
265 // di
266 dx3 = mk3d*(b1+b2) + mk23d*e*mJ13d/(mJ13d+r)*mJ63d/(mJ63d+p) + mk43d*j - mk34d*i*j/(1+j) - mphi_i*i;
267 // dj
268 dx4 = mk34d*i*j/(1+j) - (mk43d+mphi_j)*j;
269 // dp
270 dx5 = mk16d*r*j - mk61d*p - mphi_p*p;
271
272 double factor = mPhiE2F1*60.0; // convert non-dimensional d/dt s to d/dt in hours
273
274 rDY[0] = dx1*factor;
275 rDY[1] = dx2*factor;
276 rDY[2] = dx3*factor;
277 rDY[3] = dx4*factor;
278 rDY[4] = dx5*factor;
279 rDY[5] = dx6*factor;
280 rDY[6] = dx7*factor; // beta-cat allele 1
281 rDY[7] = dx8*factor; // beta-cat allele 2
282 rDY[8] = 0.0; // do not change the Wnt level
283}
284
285const boost::shared_ptr<AbstractCellMutationState> WntCellCycleOdeSystem::GetMutationState() const
286{
287 return mpMutationState;
288}
289
290bool WntCellCycleOdeSystem::CalculateStoppingEvent(double time, const std::vector<double>& rY)
291{
292 double r = rY[0];
293 double e = rY[1];
294 double p = rY[4];
295 double dY1 = mkpd+mk2d*(mad*mad+e*e)/(1+e*e)*mJ12d/(mJ12d+r)*mJ62d/(mJ62d+p) - e;
296 double factor = mPhiE2F1*60.0; // Convert non-dimensional d/dt s to d/dt in hours.
297 dY1 = dY1*factor;
298
299 assert(!std::isnan(rY[1]));
300 assert(!std::isnan(dY1));
301 return (rY[1] > 1.0 && dY1 > 0.0);
302}
303
304double WntCellCycleOdeSystem::CalculateRootFunction(double time, const std::vector<double>& rY)
305{
306 return rY[1] - 1.0;
307}
308
309template<>
311{
312 this->mVariableNames.push_back("pRb");
313 this->mVariableUnits.push_back("non_dim");
314 this->mInitialConditions.push_back(7.357000000000000e-01);
315
316 this->mVariableNames.push_back("E2F1");
317 this->mVariableUnits.push_back("non_dim");
318 this->mInitialConditions.push_back(1.713000000000000e-01);
319
320 this->mVariableNames.push_back("CycD_i");
321 this->mVariableUnits.push_back("non_dim");
322 this->mInitialConditions.push_back(6.900000000000001e-02);
323
324 this->mVariableNames.push_back("CycD_a");
325 this->mVariableUnits.push_back("non_dim");
326 this->mInitialConditions.push_back(3.333333333333334e-03);
327
328 this->mVariableNames.push_back("pRb_p");
329 this->mVariableUnits.push_back("non_dim");
330 this->mInitialConditions.push_back(1.000000000000000e-04);
331
332 this->mVariableNames.push_back("APC");
333 this->mVariableUnits.push_back("non_dim");
334 this->mInitialConditions.push_back(NAN); // will be filled in later
335
336 this->mVariableNames.push_back("Beta_Cat1");
337 this->mVariableUnits.push_back("non_dim");
338 this->mInitialConditions.push_back(NAN); // will be filled in later
339
340 this->mVariableNames.push_back("Beta_Cat2");
341 this->mVariableUnits.push_back("non_dim");
342 this->mInitialConditions.push_back(NAN); // will be filled in later
343
344 this->mVariableNames.push_back("Wnt");
345 this->mVariableUnits.push_back("non_dim");
346 this->mInitialConditions.push_back(NAN); // will be filled in later
347
348 this->mInitialised = true;
349}
350
352{
353 return mWntLevel;
354}
355
356// Serialization for Boost >= 1.36
#define CHASTE_CLASS_EXPORT(T)
void SetStateVariables(const std::vector< double > &rStateVariables)
void SetDefaultInitialCondition(unsigned index, double initialCondition)
boost::shared_ptr< AbstractOdeSystemInformation > mpSystemInfo
boost::shared_ptr< AbstractCellMutationState > mpMutationState
void SetMutationState(boost::shared_ptr< AbstractCellMutationState > pMutationState)
const boost::shared_ptr< AbstractCellMutationState > GetMutationState() const
WntCellCycleOdeSystem(double wntLevel=0.0, boost::shared_ptr< AbstractCellMutationState > pMutationState=boost::shared_ptr< AbstractCellMutationState >(), std::vector< double > stateVariables=std::vector< double >())
void EvaluateYDerivatives(double time, const std::vector< double > &rY, std::vector< double > &rDY)
double CalculateRootFunction(double time, const std::vector< double > &rY)
bool CalculateStoppingEvent(double time, const std::vector< double > &rY)