Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
VanLeeuwen2009WntSwatCellCycleOdeSystem.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 "VanLeeuwen2009WntSwatCellCycleOdeSystem.hpp"
37#include "CellwiseOdeSystemInformation.hpp"
38#include "ApcOneHitCellMutationState.hpp"
39#include "ApcTwoHitCellMutationState.hpp"
40#include "BetaCateninOneHitCellMutationState.hpp"
41
43 double wntLevel,
44 boost::shared_ptr<AbstractCellMutationState> pMutationState,
45 std::vector<double> stateVariables)
47 mpMutationState(pMutationState),
48 mHypothesis(hypothesis),
49 mWntLevel(wntLevel)
50{
51 if (hypothesis!=1 && hypothesis!=2)
52 {
53 EXCEPTION("You must set up this cell cycle ODE system with hypothesis one or two.");
54 }
55
57
84 Init(); // set up parameter values
85
86 double d_d_hat = mDd + mXiD*wntLevel;
87 double d_d_x_hat = mDdx + mXiDx*wntLevel;
88 double d_x_hat = mDx + mXiX*wntLevel;
89 double p_c_hat = mPc + mXiC*wntLevel;
90
91 double sigma_D = 0.0; // for healthy cells
92 //double sigma_B = 0.0; // for healthy cells - never used!
93
94 if (!mpMutationState)
95 {
96 // No mutations specified
97 }
99 {
100 sigma_D = 0.5;
101 }
103 {
104 sigma_D = 1.0;
105 }
107 {
108 //sigma_B = 0.5; // never used!
109 }
110 // Other mutations have no effect.
111
112 // Cell-specific initial conditions
113 double steady_D = ((1.0-sigma_D)*mSd*mSx)/((1.0-sigma_D)*mSd*d_d_hat + d_x_hat*(d_d_hat + d_d_x_hat));
114 SetDefaultInitialCondition(5, steady_D); // Destruction complex (APC/Axin/GSK3B)
115
116 double temp = (mSx*(d_d_hat+d_d_x_hat))/((1.0-sigma_D)*mSd*d_d_hat+d_x_hat*(d_d_hat+d_d_x_hat));
117 SetDefaultInitialCondition(6, temp); // Axin
118
119 double steady_Cf = ((mSc-mDc*mKd - mPu*steady_D)+sqrt(SmallPow((mSc-mDc*mKd - mPu*steady_D),2) + (4.0*mSc*mDc*mKd)))/(2.0*mDc);
120 temp = (mPu*steady_D*steady_Cf)/(mDu*(steady_Cf+mKd));
121 SetDefaultInitialCondition(7, temp); // beta-catenin to be ubiquitinated
122
123 double theta = mDc + (mPu*steady_D)/(steady_Cf + mKd);
124
125 double steady_Co = ( mSc - p_c_hat - theta*mKc + sqrt(4.0*mSc*theta*mKc + SmallPow((mSc - p_c_hat - theta*mKc),2)) )/(2.0*theta);
126 SetDefaultInitialCondition(8, steady_Co); // Open form beta-catenin
127
128 double steady_Cc = steady_Cf - steady_Co;
129
130 if ((steady_Cc < 0) && (steady_Cc+100*DBL_EPSILON > 0) ) // stop protein values going -ve
131 {
132 steady_Cc = 0.0;
133 }
134 SetDefaultInitialCondition(9, steady_Cc); // Closed form beta-catenin
135
136 SetDefaultInitialCondition(12, mSa/mDa); // 'Free' adhesion molecules
137
138 SetDefaultInitialCondition(13, mSa*mSca*steady_Co/(mDa*mDca)); // Co-A Adhesion complex
139
140 SetDefaultInitialCondition(15, mSt/mDt); // `Free' transcription molecules (TCF)
141
142 SetDefaultInitialCondition(16, mSct*mSt*steady_Co/(mDt*mDct)); // Co-T open form beta-catenin/TCF
143
144 SetDefaultInitialCondition(17, mSct*mSt*steady_Cc/(mDt*mDct)); // Cc-T closed beta-catenin/TCF
145
146 temp = (mSct*mSt*mSy*steady_Cf)/(mDy*(mSct*mSt*steady_Cf + mDct*mDt*mKt));
147 SetDefaultInitialCondition(20, temp); // Wnt target protein
148
149 SetDefaultInitialCondition(21, wntLevel); // Wnt stimulus
150
151 if (stateVariables != std::vector<double>())
152 {
153 SetStateVariables(stateVariables);
154 }
155}
156
157void VanLeeuwen2009WntSwatCellCycleOdeSystem::SetMutationState(boost::shared_ptr<AbstractCellMutationState> pMutationState)
158{
159 mpMutationState = pMutationState;
160}
161
166
168{
169 // Swat (2004) parameters
170 double k1 = 1.0;
171 double k2 = 1.6;
172 double k3 = 0.05;
173 double k16 = 0.4;
174 double k34 = 0.04;
175 double k43 = 0.01;
176 double k61 = 0.3;
177 double k23 = 0.3;
178 double a = 0.04;
179 double J11 = 0.5;
180 double J12 = 5.0;
181 double J61 = 5.0;
182 double J62 = 8.0;
183 double J13 = 0.002;
184 double J63 = 2.0;
185 double Km1 = 0.5;
186 double Km2 = 4.0;
187 double Km4 = 0.3;
188 double kp = 0.05;
189 double phi_pRb = 0.005;
190 double phi_E2F1 = 0.1;
191 double phi_CycDi = 0.023;
192 double phi_CycDa = 0.03;
193 double phi_pRbp = 0.06;
194
195 // Value of the mitogenic factor to make the van Leeuwen model influence cell cycle just the same
196 double mitogenic_factorF = 1.0/25.0;
197
198 // Non-dimensionalise parameters
199 mk2d = k2/(Km2*phi_E2F1);
200 mk3d = k3*mitogenic_factorF/(Km4*phi_E2F1);
201 mk34d = k34/phi_E2F1;
202 mk43d = k43/phi_E2F1;
203 mk23d = k23*Km2/(Km4*phi_E2F1);
204 mad = a/Km2;
205 mJ11d = J11*phi_E2F1/k1;
206 mJ12d = J12*phi_E2F1/k1;
207 mJ13d = J13*phi_E2F1/k1;
208 mJ61d = J61*phi_E2F1/k1;
209 mJ62d = J62*phi_E2F1/k1;
210 mJ63d = J63*phi_E2F1/k1;
211 mKm1d = Km1/Km2;
212 mkpd = kp/(Km2*phi_E2F1);
213 mphi_r = phi_pRb/phi_E2F1;
214 mphi_i = phi_CycDi/phi_E2F1;
215 mphi_j = phi_CycDa/phi_E2F1;
216 mphi_p = phi_pRbp/phi_E2F1;
217 mk16d = k16*Km4/phi_E2F1;
218 mk61d = k61/phi_E2F1;
219 mPhiE2F1 = phi_E2F1;
220
221 // Initialize van Leeuwen model parameters
222 mSa = 20; // nM/h
223 mSca = 250; // (nMh)^-1
224 mSc = 25; // nM/h
225 mSct = 30; // (nMh)^-1
226 mSd = 100; // h^-1
227 mSt = 10; // nM/h
228 mSx = 10; // nM/h
229 mSy = 10; // h^-1
230 mDa = 2; // h^-1
231 mDca = 350; // h^-1
232 mDc = 1; // h^-1
233 mDct = 750; // h^-1
234 mDd = 5; // h^-1
235 mDdx = 5; // h^-1
236 mDt = 0.4; // h^-1
237 mDu = 50; // h^-1
238 mDx = 100; // h^-1
239 mDy = 1; // h^-1
240 mKc = 200; // nM
241 mKd = 5; // nM
242 mKt = 50; // nM
243 mPc = 0.0; // h^-1
244 mPu = 100; // h^-1
245 mXiD = 5; // h^-1
246 mXiDx = 5; // h^-1
247 mXiX = 200; // h^-1
248 if (mHypothesis==1)
249 {
250 mXiC = 0.0; // h^-1 (FOR HYPOTHESIS ONE)
251 }
252 else
253 {
254 mXiC = 5000.0; // h^-1 (FOR HYPOTHESIS TWO)
255 }
256}
257
258void VanLeeuwen2009WntSwatCellCycleOdeSystem::EvaluateYDerivatives(double time, const std::vector<double>& rY, std::vector<double>& rDY)
259{
260 double r = rY[0];
261 double e = rY[1];
262 double i = rY[2];
263 double j = rY[3];
264 double p = rY[4];
265
266 double dx1 = 0.0;
267 double dx2 = 0.0;
268 double dx3 = 0.0;
269 double dx4 = 0.0;
270 double dx5 = 0.0;
271
272 // Bit back-to-front, but work out the Wnt section first...
273
274 // Variables
275 double D = rY[5];
276 double X = rY[6];
277 double Cu = rY[7];
278 double Co = rY[8];
279 double Cc = rY[9];
280 double Mo = rY[10];
281 double Mc = rY[11];
282 double A = rY[12];
283 double Ca = rY[13];
284 double Ma = rY[14];
285 double T = rY[15];
286 double Cot = rY[16];
287 double Cct = rY[17];
288 double Mot = rY[18];
289 double Mct = rY[19];
290 double Y = rY[20];
291 double stimulus_wnt = rY[21];
292
293 // Totals
294 double Cf = Cc+Co;
295 double Ct = Cct+Cot;
296 double Mf = Mc+Mo;
297 double Mt = Mct+Mot;
298
299 double d_d_hat = mDd + mXiD*stimulus_wnt;
300 double d_d_x_hat = mDdx + mXiDx*stimulus_wnt;
301 double d_x_hat = mDx + mXiX*stimulus_wnt;
302 double p_c_hat = mPc + mXiC*stimulus_wnt;
303
304 double sigma_D = 0.0; // for healthy cells
305 double sigma_B = 0.0; // for healthy cells
306
307 if (!mpMutationState)
308 {
309 // No mutations specified
310 }
312 {
313 sigma_D = 0.5;
314 }
316 {
317 sigma_D = 1.0;
318 }
320 {
321 sigma_B = 0.5;
322 }
323 // Other mutations have no effect.
324
325 // Now the cell cycle stuff...
326
327 // dr
328 dx1 = e/(mKm1d+e)*mJ11d/(mJ11d+r)*mJ61d/(mJ61d+p) - mk16d*r*j+mk61d*p-mphi_r*r;
329 // de
330 dx2 = mkpd+mk2d*(mad*mad+e*e)/(1+e*e)*mJ12d/(mJ12d+r)*mJ62d/(mJ62d+p) - e;
331 // di - changed to include Ct+Mt - transcriptional beta-catenin
332 dx3 = mk3d*(Ct+Mt) + mk23d*e*mJ13d/(mJ13d+r)*mJ63d/(mJ63d+p) + mk43d*j - mk34d*i*j/(1+j) - mphi_i*i;
333 // dj
334 dx4 = mk34d*i*j/(1+j) - (mk43d+mphi_j)*j;
335 // dp
336 dx5 = mk16d*r*j - mk61d*p - mphi_p*p;
337
338 double factor = mPhiE2F1*60.0; // Convert non-dimensional d/dt s to d/dt in hours.
339
340 rDY[0] = dx1*factor;
341 rDY[1] = dx2*factor;
342 rDY[2] = dx3*factor;
343 rDY[3] = dx4*factor;
344 rDY[4] = dx5*factor;
345
346 // The van Leeuwen ODE system
347 rDY[5] = (1.0-sigma_D)*mSd*X - (d_d_hat + d_d_x_hat)*D;
348 rDY[6] = mSx - (1.0-sigma_D)*mSd*X - d_x_hat*X + d_d_x_hat*D;
349 rDY[7] = (mPu*D*Cf)/(Cf+mKd) - mDu*Cu;
350
351 rDY[8] = (1.0-sigma_B)*mSc + mDca*Ca + mDct*Cot - (mSca*A + mSct*T + mDc)*Co
352 - (p_c_hat*Co)/(Co + Mo + mKc) - (mPu*D*Co)/(Cf+mKd);
353
354 rDY[9] = (p_c_hat*Co)/(Co + Mo + mKc) + mDct*Cct - (mSct*T + mDc)*Cc
355 - (mPu*D*Cc)/(Cf+mKd);
356
357 rDY[10] = sigma_B*mSc + mDca*Ma + mDct*Mot - (mSca*A + mSct*T + mDc)*Mo
358 - (p_c_hat*Mo)/(Co + Mo + mKc);
359
360 rDY[11] = (p_c_hat*Mo)/(Co + Mo + mKc) + mDct*Mct - (mSct*T + mDc)*Mc;
361 rDY[12] = mSa + mDca*(Ca+Ma) - (mSca*(Co+Mo) + mDa)*A;
362 rDY[13] = mSca*Co*A - mDca*Ca;
363 rDY[14] = mSca*Mo*A - mDca*Ma;
364 rDY[15] = mSt + mDct*(Ct+Mt) - mSct*(Cf+Mf)*T - mDt*T;
365 rDY[16] = mSct*Co*T - mDct*Cot;
366 rDY[17] = mSct*Cc*T - mDct*Cct;
367 rDY[18] = mSct*Mo*T - mDct*Mot;
368 rDY[19] = mSct*Mc*T - mDct*Mct;
369 rDY[20] = (mSy*(Ct+Mt))/(Ct + Mt + mKt) - mDy*Y;
370 rDY[21] = 0.0; // don't interfere with Wnt stimulus
371}
372
373const boost::shared_ptr<AbstractCellMutationState> VanLeeuwen2009WntSwatCellCycleOdeSystem::GetMutationState() const
374{
375 return mpMutationState;
376}
377
378bool VanLeeuwen2009WntSwatCellCycleOdeSystem::CalculateStoppingEvent(double time, const std::vector<double>& rY)
379{
380 std::vector<double> dy(rY.size());
381 EvaluateYDerivatives(time, rY, dy);
382
383 return (rY[1] > 1.0 && dy[1] > 0.0);
384}
385
386double VanLeeuwen2009WntSwatCellCycleOdeSystem::CalculateRootFunction(double time, const std::vector<double>& rY)
387{
388 return rY[1] - 1.0;
389}
390
391template<>
393{
394 this->mVariableNames.push_back("pRb");
395 this->mVariableUnits.push_back("non_dim");
396 this->mInitialConditions.push_back(7.357000000000000e-01);
397
398 this->mVariableNames.push_back("E2F1");
399 this->mVariableUnits.push_back("non_dim");
400 this->mInitialConditions.push_back(1.713000000000000e-01);
401
402 this->mVariableNames.push_back("CycD_i");
403 this->mVariableUnits.push_back("non_dim");
404 this->mInitialConditions.push_back(6.900000000000001e-02);
405
406 this->mVariableNames.push_back("CycD_a");
407 this->mVariableUnits.push_back("non_dim");
408 this->mInitialConditions.push_back(3.333333333333334e-03);
409
410 this->mVariableNames.push_back("pRb_p");
411 this->mVariableUnits.push_back("non_dim");
412 this->mInitialConditions.push_back(1.000000000000000e-04);
413
414 this->mVariableNames.push_back("D"); // Destruction complex (APC/Axin/GSK3B)
415 this->mVariableUnits.push_back("nM");
416 this->mInitialConditions.push_back(NAN); // will be filled in later
417
418 this->mVariableNames.push_back("X"); // Axin
419 this->mVariableUnits.push_back("nM");
420 this->mInitialConditions.push_back(NAN); // will be filled in later
421
422 this->mVariableNames.push_back("Cu"); // beta-catenin to be ubiquitinated
423 this->mVariableUnits.push_back("nM");
424 this->mInitialConditions.push_back(NAN); // will be filled in later
425
426 this->mVariableNames.push_back("Co"); // Open form beta-catenin
427 this->mVariableUnits.push_back("nM");
428 this->mInitialConditions.push_back(NAN); // will be filled in later
429
430 this->mVariableNames.push_back("Cc"); // Closed form beta-catenin
431 this->mVariableUnits.push_back("nM");
432 this->mInitialConditions.push_back(NAN); // will be filled in later
433
434 this->mVariableNames.push_back("Mo"); // Open form mutant beta-catenin
435 this->mVariableUnits.push_back("nM");
436 this->mInitialConditions.push_back(0.0);
437
438 this->mVariableNames.push_back("Mc"); // Closed form mutant beta-catenin
439 this->mVariableUnits.push_back("nM");
440 this->mInitialConditions.push_back(0.0);
441
442 this->mVariableNames.push_back("A"); // 'Free' adhesion molecules
443 this->mVariableUnits.push_back("nM");
444 this->mInitialConditions.push_back(NAN); // will be filled in later
445
446 this->mVariableNames.push_back("Ca"); // Co-A Adhesion complex
447 this->mVariableUnits.push_back("nM");
448 this->mInitialConditions.push_back(NAN); // will be filled in later
449
450 this->mVariableNames.push_back("Ma"); // Mo-A Mutant adhesion complex
451 this->mVariableUnits.push_back("nM");
452 this->mInitialConditions.push_back(0.0);
453
454 this->mVariableNames.push_back("T"); // `Free' transcription molecules (TCF)
455 this->mVariableUnits.push_back("nM");
456 this->mInitialConditions.push_back(NAN); // will be filled in later
457
458 this->mVariableNames.push_back("Cot"); // Co-T open form beta-catenin/TCF
459 this->mVariableUnits.push_back("nM");
460 this->mInitialConditions.push_back(NAN); // will be filled in later
461
462 this->mVariableNames.push_back("Cct"); // Cc-T closed beta-catenin/TCF
463 this->mVariableUnits.push_back("nM");
464 this->mInitialConditions.push_back(NAN); // will be filled in later
465
466 this->mVariableNames.push_back("Mot"); // Mo-T open form mutant beta-catenin/TCF
467 this->mVariableUnits.push_back("nM");
468 this->mInitialConditions.push_back(0.0);
469
470 this->mVariableNames.push_back("Mct"); // Mc-T closed form mutant beta-catenin/TCF
471 this->mVariableUnits.push_back("nM");
472 this->mInitialConditions.push_back(0.0);
473
474 this->mVariableNames.push_back("Y"); // Wnt target protein
475 this->mVariableUnits.push_back("nM");
476 this->mInitialConditions.push_back(NAN); // will be filled in later
477
478 this->mVariableNames.push_back("Sw"); // Wnt stimulus
479 this->mVariableUnits.push_back("nM");
480 this->mInitialConditions.push_back(NAN); // will be filled in later
481
482 this->mInitialised = true;
483}
484
489
494
495// Serialization for Boost >= 1.36
#define EXCEPTION(message)
#define CHASTE_CLASS_EXPORT(T)
void SetStateVariables(const std::vector< double > &rStateVariables)
void SetDefaultInitialCondition(unsigned index, double initialCondition)
boost::shared_ptr< AbstractOdeSystemInformation > mpSystemInfo
const boost::shared_ptr< AbstractCellMutationState > GetMutationState() const
double CalculateRootFunction(double time, const std::vector< double > &rY)
bool CalculateStoppingEvent(double time, const std::vector< double > &rY)
void EvaluateYDerivatives(double time, const std::vector< double > &rY, std::vector< double > &rDY)
void SetMutationState(boost::shared_ptr< AbstractCellMutationState > pMutationState)
boost::shared_ptr< AbstractCellMutationState > mpMutationState
VanLeeuwen2009WntSwatCellCycleOdeSystem(unsigned hypothesis, double wntLevel=0.0, boost::shared_ptr< AbstractCellMutationState > pMutationState=boost::shared_ptr< AbstractCellMutationState >(), std::vector< double > stateVariables=std::vector< double >())