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