VanLeeuwen2009WntSwatCellCycleOdeSystem.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 
00029 #include "VanLeeuwen2009WntSwatCellCycleOdeSystem.hpp"
00030 #include "CellwiseOdeSystemInformation.hpp"
00031 
00032 VanLeeuwen2009WntSwatCellCycleOdeSystem::VanLeeuwen2009WntSwatCellCycleOdeSystem(unsigned hypothesis,
00033                                                                                  double wntLevel,
00034                                                                                  boost::shared_ptr<AbstractCellMutationState> pMutationState,
00035                                                                                  std::vector<double> stateVariables)
00036     : AbstractOdeSystem(22),
00037       mpMutationState(pMutationState),
00038       mHypothesis(hypothesis),
00039       mWntLevel(wntLevel)
00040 {
00041     if (hypothesis!=1 && hypothesis!=2)
00042     {
00043         EXCEPTION("You must set up this cell cycle ODE system with hypothesis one or two.");
00044     }
00045 
00046     mpSystemInfo.reset(new CellwiseOdeSystemInformation<VanLeeuwen2009WntSwatCellCycleOdeSystem>);
00047 
00074     Init(); // set up parameter values
00075 
00076     double d_d_hat = mDd + mXiD*wntLevel;
00077     double d_d_x_hat = mDdx + mXiDx*wntLevel;
00078     double d_x_hat = mDx + mXiX*wntLevel;
00079     double p_c_hat = mPc + mXiC*wntLevel;
00080 
00081     double sigma_D = 0.0; // for healthy cells
00082     double sigma_B = 0.0; // for healthy cells
00083 
00084     if (!mpMutationState)
00085     {
00086         // No mutations specified
00087     }
00088     else if (mpMutationState->IsType<ApcOneHitCellMutationState>())
00089     {
00090         sigma_D = 0.5;
00091     }
00092     else if (mpMutationState->IsType<ApcTwoHitCellMutationState>())
00093     {
00094         sigma_D = 1.0;
00095     }
00096     else if (mpMutationState->IsType<BetaCateninOneHitCellMutationState>())
00097     {
00098         sigma_B = 0.5;
00099     }
00100     // Other mutations have no effect.
00101 
00102     // Cell-specific initial conditions
00103     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));
00104     SetDefaultInitialCondition(5, steady_D); // Destruction complex (APC/Axin/GSK3B)
00105 
00106     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));
00107     SetDefaultInitialCondition(6, temp);  // Axin
00108 
00109     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);
00110     temp = (mPu*steady_D*steady_Cf)/(mDu*(steady_Cf+mKd));
00111     SetDefaultInitialCondition(7, temp); // beta-catenin to be ubiquitinated
00112 
00113     double theta = mDc + (mPu*steady_D)/(steady_Cf + mKd);
00114 
00115     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);
00116     SetDefaultInitialCondition(8, steady_Co); // Open form beta-catenin
00117 
00118     double steady_Cc = steady_Cf - steady_Co;
00119 
00120     if ((steady_Cc < 0) && (steady_Cc+100*DBL_EPSILON > 0) ) // stop protein values going -ve
00121     {
00122         steady_Cc = 0.0;
00123     }
00124     SetDefaultInitialCondition(9, steady_Cc); // Closed form beta-catenin
00125 
00126     SetDefaultInitialCondition(12, mSa/mDa);  // 'Free' adhesion molecules
00127 
00128     SetDefaultInitialCondition(13, mSa*mSca*steady_Co/(mDa*mDca)); // Co-A Adhesion complex
00129 
00130     SetDefaultInitialCondition(15, mSt/mDt); // `Free' transcription molecules (TCF)
00131 
00132     SetDefaultInitialCondition(16, mSct*mSt*steady_Co/(mDt*mDct)); // Co-T open form beta-catenin/TCF
00133 
00134     SetDefaultInitialCondition(17, mSct*mSt*steady_Cc/(mDt*mDct)); // Cc-T closed beta-catenin/TCF
00135 
00136     temp = (mSct*mSt*mSy*steady_Cf)/(mDy*(mSct*mSt*steady_Cf + mDct*mDt*mKt));
00137     SetDefaultInitialCondition(20, temp); // Wnt target protein
00138 
00139     SetDefaultInitialCondition(21, wntLevel); // Wnt stimulus
00140 
00141     if (stateVariables != std::vector<double>())
00142     {
00143         SetStateVariables(stateVariables);
00144     }
00145 }
00146 
00147 void VanLeeuwen2009WntSwatCellCycleOdeSystem::SetMutationState(boost::shared_ptr<AbstractCellMutationState> pMutationState)
00148 {
00149     mpMutationState = pMutationState;
00150 }
00151 
00152 VanLeeuwen2009WntSwatCellCycleOdeSystem::~VanLeeuwen2009WntSwatCellCycleOdeSystem()
00153 {
00154     // Do nothing
00155 }
00156 
00157 void VanLeeuwen2009WntSwatCellCycleOdeSystem::Init()
00158 {
00159     // Swat (2004) parameters
00160     double k1 = 1.0;
00161     double k2 = 1.6;
00162     double k3 = 0.05;
00163     double k16 = 0.4;
00164     double k34 = 0.04;
00165     double k43 = 0.01;
00166     double k61 = 0.3;
00167     double k23 = 0.3;
00168     double a = 0.04;
00169     double J11 = 0.5;
00170     double J12 = 5.0;
00171     double J61 = 5.0;
00172     double J62 = 8.0;
00173     double J13 = 0.002;
00174     double J63 = 2.0;
00175     double Km1 = 0.5;
00176     double Km2 = 4.0;
00177     double Km4 = 0.3;
00178     double kp = 0.05;
00179     double phi_pRb = 0.005;
00180     double phi_E2F1 = 0.1;
00181     double phi_CycDi = 0.023;
00182     double phi_CycDa = 0.03;
00183     double phi_pRbp = 0.06;
00184 
00185     // Value of the mitogenic factor to make the van Leeuwen model influence cell cycle just the same
00186     double mitogenic_factorF = 1.0/25.0;
00187 
00188     // Non-dimensionalise parameters
00189     mk2d = k2/(Km2*phi_E2F1);
00190     mk3d = k3*mitogenic_factorF/(Km4*phi_E2F1);
00191     mk34d = k34/phi_E2F1;
00192     mk43d = k43/phi_E2F1;
00193     mk23d = k23*Km2/(Km4*phi_E2F1);
00194     mad = a/Km2;
00195     mJ11d = J11*phi_E2F1/k1;
00196     mJ12d = J12*phi_E2F1/k1;
00197     mJ13d = J13*phi_E2F1/k1;
00198     mJ61d = J61*phi_E2F1/k1;
00199     mJ62d = J62*phi_E2F1/k1;
00200     mJ63d = J63*phi_E2F1/k1;
00201     mKm1d = Km1/Km2;
00202     mkpd = kp/(Km2*phi_E2F1);
00203     mphi_r = phi_pRb/phi_E2F1;
00204     mphi_i = phi_CycDi/phi_E2F1;
00205     mphi_j = phi_CycDa/phi_E2F1;
00206     mphi_p = phi_pRbp/phi_E2F1;
00207     mk16d = k16*Km4/phi_E2F1;
00208     mk61d = k61/phi_E2F1;
00209     mPhiE2F1 = phi_E2F1;
00210 
00211     // Initialize van Leeuwen model parameters
00212     mSa = 20;   //  nM/h
00213     mSca = 250; //  (nMh)^-1
00214     mSc = 25;   //  nM/h
00215     mSct = 30;  //  (nMh)^-1
00216     mSd = 100;  //  h^-1
00217     mSt = 10;   //  nM/h
00218     mSx = 10;   //  nM/h
00219     mSy = 10;   //  h^-1
00220     mDa = 2;    //  h^-1
00221     mDca = 350; //  h^-1
00222     mDc = 1;    //  h^-1
00223     mDct = 750; //  h^-1
00224     mDd = 5;    //  h^-1
00225     mDdx = 5;   //  h^-1
00226     mDt = 0.4;  //  h^-1
00227     mDu = 50;   //  h^-1
00228     mDx = 100;  //  h^-1
00229     mDy = 1;    //  h^-1
00230     mKc = 200;  //  nM
00231     mKd = 5;    //  nM
00232     mKt = 50;   //  nM
00233     mPc = 0.0;  //  h^-1
00234     mPu = 100;  //  h^-1
00235     mXiD = 5;   //  h^-1
00236     mXiDx = 5;  //  h^-1
00237     mXiX = 200; //  h^-1
00238     if (mHypothesis==1)
00239     {
00240         mXiC = 0.0; //  h^-1 (FOR HYPOTHESIS ONE)
00241     }
00242     else
00243     {
00244         mXiC = 5000.0;   //  h^-1 (FOR HYPOTHESIS TWO)
00245     }
00246 }
00247 
00248 void VanLeeuwen2009WntSwatCellCycleOdeSystem::EvaluateYDerivatives(double time, const std::vector<double>& rY, std::vector<double>& rDY)
00249 {
00250     double r = rY[0];
00251     double e = rY[1];
00252     double i = rY[2];
00253     double j = rY[3];
00254     double p = rY[4];
00255 
00256     double dx1 = 0.0;
00257     double dx2 = 0.0;
00258     double dx3 = 0.0;
00259     double dx4 = 0.0;
00260     double dx5 = 0.0;
00261 
00262     // Bit back-to-front, but work out the Wnt section first...
00263 
00264     // Variables
00265     double D = rY[5];
00266     double X = rY[6];
00267     double Cu = rY[7];
00268     double Co = rY[8];
00269     double Cc = rY[9];
00270     double Mo = rY[10];
00271     double Mc = rY[11];
00272     double A = rY[12];
00273     double Ca = rY[13];
00274     double Ma = rY[14];
00275     double T = rY[15];
00276     double Cot = rY[16];
00277     double Cct = rY[17];
00278     double Mot = rY[18];
00279     double Mct = rY[19];
00280     double Y = rY[20];
00281     double stimulus_wnt = rY[21];
00282 
00283     // Totals
00284     double Cf = Cc+Co;
00285     double Ct = Cct+Cot;
00286     double Mf = Mc+Mo;
00287     double Mt = Mct+Mot;
00288 
00289     double d_d_hat = mDd + mXiD*stimulus_wnt;
00290     double d_d_x_hat = mDdx + mXiDx*stimulus_wnt;
00291     double d_x_hat = mDx + mXiX*stimulus_wnt;
00292     double p_c_hat = mPc + mXiC*stimulus_wnt;
00293 
00294     double sigma_D = 0.0;   // for healthy cells
00295     double sigma_B = 0.0;   // for healthy cells
00296 
00297     if (!mpMutationState)
00298     {
00299         // No mutations specified
00300     }
00301     else if (mpMutationState->IsType<ApcOneHitCellMutationState>())
00302     {
00303         sigma_D = 0.5;
00304     }
00305     else if (mpMutationState->IsType<ApcTwoHitCellMutationState>())
00306     {
00307         sigma_D = 1.0;
00308     }
00309     else if (mpMutationState->IsType<BetaCateninOneHitCellMutationState>())
00310     {
00311         sigma_B = 0.5;
00312     }
00313     // Other mutations have no effect.
00314 
00315     // Now the cell cycle stuff...
00316 
00317     // dr
00318     dx1 = e/(mKm1d+e)*mJ11d/(mJ11d+r)*mJ61d/(mJ61d+p) - mk16d*r*j+mk61d*p-mphi_r*r;
00319     // de
00320     dx2 = mkpd+mk2d*(mad*mad+e*e)/(1+e*e)*mJ12d/(mJ12d+r)*mJ62d/(mJ62d+p) - e;
00321     // di - changed to include Ct+Mt - transcriptional beta-catenin
00322     dx3 = mk3d*(Ct+Mt) + mk23d*e*mJ13d/(mJ13d+r)*mJ63d/(mJ63d+p) + mk43d*j - mk34d*i*j/(1+j) - mphi_i*i;
00323     // dj
00324     dx4 = mk34d*i*j/(1+j) - (mk43d+mphi_j)*j;
00325     // dp
00326     dx5 = mk16d*r*j - mk61d*p - mphi_p*p;
00327 
00328     double factor = mPhiE2F1*60.0;  // Convert non-dimensional d/dt s to d/dt in hours.
00329 
00330     rDY[0] = dx1*factor;
00331     rDY[1] = dx2*factor;
00332     rDY[2] = dx3*factor;
00333     rDY[3] = dx4*factor;
00334     rDY[4] = dx5*factor;
00335 
00336     // The van Leeuwen ODE system
00337     rDY[5] = (1.0-sigma_D)*mSd*X - (d_d_hat + d_d_x_hat)*D;
00338     rDY[6] = mSx - (1.0-sigma_D)*mSd*X - d_x_hat*X + d_d_x_hat*D;
00339     rDY[7] = (mPu*D*Cf)/(Cf+mKd) - mDu*Cu;
00340 
00341     rDY[8] = (1.0-sigma_B)*mSc + mDca*Ca + mDct*Cot - (mSca*A + mSct*T + mDc)*Co
00342              - (p_c_hat*Co)/(Co + Mo + mKc) - (mPu*D*Co)/(Cf+mKd);
00343 
00344     rDY[9] = (p_c_hat*Co)/(Co + Mo + mKc) + mDct*Cct - (mSct*T + mDc)*Cc
00345              - (mPu*D*Cc)/(Cf+mKd);
00346 
00347     rDY[10] = sigma_B*mSc + mDca*Ma + mDct*Mot - (mSca*A + mSct*T + mDc)*Mo
00348              - (p_c_hat*Mo)/(Co + Mo + mKc);
00349 
00350     rDY[11] = (p_c_hat*Mo)/(Co + Mo + mKc) + mDct*Mct - (mSct*T + mDc)*Mc;
00351     rDY[12] = mSa + mDca*(Ca+Ma) - (mSca*(Co+Mo) + mDa)*A;
00352     rDY[13] = mSca*Co*A - mDca*Ca;
00353     rDY[14] = mSca*Mo*A - mDca*Ma;
00354     rDY[15] = mSt + mDct*(Ct+Mt) - mSct*(Cf+Mf)*T - mDt*T;
00355     rDY[16] = mSct*Co*T - mDct*Cot;
00356     rDY[17] = mSct*Cc*T - mDct*Cct;
00357     rDY[18] = mSct*Mo*T - mDct*Mot;
00358     rDY[19] = mSct*Mc*T - mDct*Mct;
00359     rDY[20] = (mSy*(Ct+Mt))/(Ct + Mt + mKt) - mDy*Y;
00360     rDY[21] = 0.0;  // don't interfere with Wnt stimulus
00361 }
00362 
00363 const boost::shared_ptr<AbstractCellMutationState> VanLeeuwen2009WntSwatCellCycleOdeSystem::GetMutationState() const
00364 {
00365     return mpMutationState;
00366 }
00367 
00368 bool VanLeeuwen2009WntSwatCellCycleOdeSystem::CalculateStoppingEvent(double time, const std::vector<double>& rY)
00369 {
00370     std::vector<double> dy(rY.size());
00371     EvaluateYDerivatives(time, rY, dy);
00372 
00373     return (rY[1] > 1.0 && dy[1] > 0.0);
00374 }
00375 
00376 double VanLeeuwen2009WntSwatCellCycleOdeSystem::CalculateRootFunction(double time, const std::vector<double>& rY)
00377 {
00378     return rY[1] - 1.0;
00379 }
00380 
00381 template<>
00382 void CellwiseOdeSystemInformation<VanLeeuwen2009WntSwatCellCycleOdeSystem>::Initialise()
00383 {
00384     this->mVariableNames.push_back("pRb");
00385     this->mVariableUnits.push_back("non_dim");
00386     this->mInitialConditions.push_back(7.357000000000000e-01);
00387 
00388     this->mVariableNames.push_back("E2F1");
00389     this->mVariableUnits.push_back("non_dim");
00390     this->mInitialConditions.push_back(1.713000000000000e-01);
00391 
00392     this->mVariableNames.push_back("CycD_i");
00393     this->mVariableUnits.push_back("non_dim");
00394     this->mInitialConditions.push_back(6.900000000000001e-02);
00395 
00396     this->mVariableNames.push_back("CycD_a");
00397     this->mVariableUnits.push_back("non_dim");
00398     this->mInitialConditions.push_back(3.333333333333334e-03);
00399 
00400     this->mVariableNames.push_back("pRb_p");
00401     this->mVariableUnits.push_back("non_dim");
00402     this->mInitialConditions.push_back(1.000000000000000e-04);
00403 
00404     this->mVariableNames.push_back("D");  // Destruction complex (APC/Axin/GSK3B)
00405     this->mVariableUnits.push_back("nM");
00406     this->mInitialConditions.push_back(NAN); // will be filled in later
00407 
00408     this->mVariableNames.push_back("X");  // Axin
00409     this->mVariableUnits.push_back("nM");
00410     this->mInitialConditions.push_back(NAN); // will be filled in later
00411 
00412     this->mVariableNames.push_back("Cu"); // beta-catenin to be ubiquitinated
00413     this->mVariableUnits.push_back("nM");
00414     this->mInitialConditions.push_back(NAN); // will be filled in later
00415 
00416     this->mVariableNames.push_back("Co"); // Open form beta-catenin
00417     this->mVariableUnits.push_back("nM");
00418     this->mInitialConditions.push_back(NAN); // will be filled in later
00419 
00420     this->mVariableNames.push_back("Cc"); // Closed form beta-catenin
00421     this->mVariableUnits.push_back("nM");
00422     this->mInitialConditions.push_back(NAN); // will be filled in later
00423 
00424     this->mVariableNames.push_back("Mo"); // Open form mutant beta-catenin
00425     this->mVariableUnits.push_back("nM");
00426     this->mInitialConditions.push_back(0.0);
00427 
00428     this->mVariableNames.push_back("Mc"); // Closed form mutant beta-catenin
00429     this->mVariableUnits.push_back("nM");
00430     this->mInitialConditions.push_back(0.0);
00431 
00432     this->mVariableNames.push_back("A");  // 'Free' adhesion molecules
00433     this->mVariableUnits.push_back("nM");
00434     this->mInitialConditions.push_back(NAN); // will be filled in later
00435 
00436     this->mVariableNames.push_back("Ca"); // Co-A Adhesion complex
00437     this->mVariableUnits.push_back("nM");
00438     this->mInitialConditions.push_back(NAN); // will be filled in later
00439 
00440     this->mVariableNames.push_back("Ma"); // Mo-A Mutant adhesion complex
00441     this->mVariableUnits.push_back("nM");
00442     this->mInitialConditions.push_back(0.0);
00443 
00444     this->mVariableNames.push_back("T"); // `Free' transcription molecules (TCF)
00445     this->mVariableUnits.push_back("nM");
00446     this->mInitialConditions.push_back(NAN); // will be filled in later
00447 
00448     this->mVariableNames.push_back("Cot"); // Co-T open form beta-catenin/TCF
00449     this->mVariableUnits.push_back("nM");
00450     this->mInitialConditions.push_back(NAN); // will be filled in later
00451 
00452     this->mVariableNames.push_back("Cct"); // Cc-T closed beta-catenin/TCF
00453     this->mVariableUnits.push_back("nM");
00454     this->mInitialConditions.push_back(NAN); // will be filled in later
00455 
00456     this->mVariableNames.push_back("Mot"); // Mo-T open form mutant beta-catenin/TCF
00457     this->mVariableUnits.push_back("nM");
00458     this->mInitialConditions.push_back(0.0);
00459 
00460     this->mVariableNames.push_back("Mct"); // Mc-T closed form mutant beta-catenin/TCF
00461     this->mVariableUnits.push_back("nM");
00462     this->mInitialConditions.push_back(0.0);
00463 
00464     this->mVariableNames.push_back("Y"); // Wnt target protein
00465     this->mVariableUnits.push_back("nM");
00466     this->mInitialConditions.push_back(NAN); // will be filled in later
00467 
00468     this->mVariableNames.push_back("Sw"); // Wnt stimulus
00469     this->mVariableUnits.push_back("nM");
00470     this->mInitialConditions.push_back(NAN); // will be filled in later
00471 
00472     this->mInitialised = true;
00473 }
00474 
00475 double VanLeeuwen2009WntSwatCellCycleOdeSystem::GetWntLevel() const
00476 {
00477     return mWntLevel;
00478 }
00479 
00480 unsigned VanLeeuwen2009WntSwatCellCycleOdeSystem::GetHypothesis() const
00481 {
00482     return mHypothesis;
00483 }
00484 
00485 // Serialization for Boost >= 1.36
00486 #include "SerializationExportWrapperForCpp.hpp"
00487 CHASTE_CLASS_EXPORT(VanLeeuwen2009WntSwatCellCycleOdeSystem)

Generated on Mon Apr 18 11:35:28 2011 for Chaste by  doxygen 1.5.5