TysonNovak2001OdeSystem.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2009
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 #include "TysonNovak2001OdeSystem.hpp"
00029 #include "OdeSystemInformation.hpp"
00030 
00031 TysonNovak2001OdeSystem::TysonNovak2001OdeSystem()
00032         : AbstractOdeSystemWithAnalyticJacobian(6)
00033 {
00034     mpSystemInfo = OdeSystemInformation<TysonNovak2001OdeSystem>::Instance();
00035     
00036     Init();
00037 }
00038 
00039 TysonNovak2001OdeSystem::~TysonNovak2001OdeSystem(void)
00040 {
00041     // Do nothing
00042 }
00043 
00044 void TysonNovak2001OdeSystem::Init()
00045 {
00046     // Initialise model parameter values
00047     mK1 = 0.04;
00048     mK2d = 0.04;
00049     mK2dd = 1.0;
00050     mK2ddd = 1.0;
00051     mCycB_threshold = 0.1;
00052     mK3d = 1.0;
00053     mK3dd = 10.0;
00054     mK4d = 2.0;
00055     mK4 = 35;
00056     mJ3 = 0.04;
00057     mJ4 = 0.04;
00058     mK5d = 0.005;
00059     mK5dd = 0.2;
00060     mK6 = 0.1;
00061     mJ5 = 0.3;
00062     mN = 4;
00063     mK7 = 1.0;
00064     mK8 = 0.5;
00065     mJ7 = 1e-3;
00066     mJ8 = 1e-3;
00067     mMad = 1.0;
00068     mK9 = 0.1;
00069     mK10 = 0.02;
00070     mK11 = 1.0;
00071     mK12d = 0.2;
00072     mK12dd = 50.0;
00073     mK12ddd = 100.0;
00074     mKeq = 1e3;
00075     mK13 = 1.0;
00076     mK14 = 1.0;
00077     mK15d = 1.5;
00078     mK15dd = 0.05;
00079     mK16d = 1.0;
00080     mK16dd = 3.0;
00081     mJ15 = 0.01;
00082     mJ16 = 0.01;
00083     mMu = 0.01;
00084     mMstar = 10.0;
00085 }
00086 
00087 void TysonNovak2001OdeSystem::EvaluateYDerivatives(double time, const std::vector<double> &rY, std::vector<double> &rDY)
00088 {
00089     double x1 = rY[0];
00090     double x2 = rY[1];
00091     double x3 = rY[2];
00092     double x4 = rY[3];
00093     double x5 = rY[4];
00094     double x6 = rY[5];
00095     
00096     double dx1 = 0.0;
00097     double dx2 = 0.0;
00098     double dx3 = 0.0;
00099     double dx4 = 0.0;
00100     double dx5 = 0.0;
00101     double dx6 = 0.0;
00102 
00103     double temp1 = 0.0;
00104     double temp2 = 0.0;
00105     double temp3 = 0.0;
00106     
00116     dx1 = mK1-(mK2d+mK2dd*x2)*x1;
00117     
00118     // The commented line below models the start transition, no cycling, without Cdc20A
00119 //    temp1 = ((mK3d)*(1.0-x2))/(mJ3+1.0-x2);
00120     
00121     temp1 = ((mK3d+mK3dd*x4)*(1.0-x2))/(mJ3+1.0-x2);
00122     temp2 = (mK4*x6*x1*x2)/(mJ4+x2);
00123     dx2 = temp1-temp2;
00124     
00125     temp1 = mK5dd*(pow(x1*x6/mJ5,mN)/(1+pow(x1*x6/mJ5,mN)));
00126     temp2 = mK6*x3;
00127     dx3 = mK5d + temp1 - temp2;
00128     
00129     temp1 = (mK7*x5*(x3-x4))/(mJ7+x3-x4);
00130     temp2 = (mK8*mMad*x4)/(mJ8+x4);
00131     temp3 = mK6*x4;
00132     dx4 = temp1 - temp2 - temp3;
00133     
00134     dx5 = mK9*x6*x1*(1.0-x5) - mK10*x5;
00135     
00136     dx6 = mMu*x6*(1.0-x6/mMstar);
00137     
00138     // Multiply by 60 beacuase the Tyson and Noval 2001 paper has time in minutes, not hours
00139     rDY[0] = dx1*60.0;
00140     rDY[1] = dx2*60.0;
00141     rDY[2] = dx3*60.0;
00142     rDY[3] = dx4*60.0;
00143     rDY[4] = dx5*60.0;
00144     rDY[5] = dx6*60.0;
00145 }
00146 
00147 void TysonNovak2001OdeSystem::AnalyticJacobian(const std::vector<double> &solutionGuess, double** jacobian, double time, double timeStep)
00148 {
00149     timeStep *= 60.0; // to scale Jacobian so in hours not minutes
00150     double x1 = solutionGuess[0];
00151     double x2 = solutionGuess[1];
00152     double x3 = solutionGuess[2];
00153     double x4 = solutionGuess[3];
00154     double x5 = solutionGuess[4];
00155     double x6 = solutionGuess[5];
00156     
00157     // f1
00158     double df1_dx1 = -mK2d - mK2dd*x2;
00159     double df1_dx2 = -mK2dd*x1;
00160     
00161     jacobian[0][0] =  1-timeStep*df1_dx1;
00162     jacobian[0][1] = -timeStep*df1_dx2;
00163     
00164     // f2
00165     double df2_dx1 = -mK4*x6*x2/(mJ4+x2);
00166     double df2_dx2 = -mJ3*(mK3d + mK3dd*x4)/(pow((mJ3 + 1 - x2),2))
00167                      -mJ4*mK4*x6*x1/(pow((mJ4+x2),2));
00168     double df2_dx4 =  mK3dd*(1-x2)/(mJ3+1-x2);
00169     double df2_dx6 = -mK4*x1*x2/(mJ4+x2);
00170     
00171     jacobian[1][0] = -timeStep*df2_dx1;
00172     jacobian[1][1] =  1-timeStep*df2_dx2;
00173     jacobian[1][3] = -timeStep*df2_dx4;
00174     jacobian[1][5] = -timeStep*df2_dx6;
00175     
00176     //f3
00177     double z = x1*x6/mJ5;
00178     double df3_dx1 = (mK5dd*x6/mJ5)*mN*pow(z,mN-1)/(pow((1-pow(z,mN)),2));
00179     double df3_dx3 = -mK6;
00180     double df3_dx6 = (mK5dd*x1/mJ5)*mN*pow(z,mN-1)/(pow((1-pow(z,mN)),2));
00181     
00182     jacobian[2][0] = -timeStep*df3_dx1;
00183     jacobian[2][2] = 1-timeStep*df3_dx3;
00184     jacobian[2][5] = -timeStep*df3_dx6;
00185     
00186     // f4
00187     double df4_dx3 =  mJ7*mK7*x5/(pow(mJ7+x3-x4,2));
00188     double df4_dx4 = -mJ7*mK7*x5/(pow(mJ7+x3-x4,2)) - mK6 - mJ8*mK8*mMad/(pow(mJ8+x4,2));
00189     double df4_dx5 =  mK7*(x3-x4)/(mJ7+x3-x4);
00190     
00191     jacobian[3][2] = -timeStep*df4_dx3;
00192     jacobian[3][3] = 1-timeStep*df4_dx4;
00193     jacobian[3][4] = -timeStep*df4_dx5;
00194     
00195     // f5
00196     double df5_dx1 =  mK9*x6*(1-x5);
00197     double df5_dx5 = -mK10 - mK9*x6*x1;
00198     double df5_dx6 =  mK9*x1*(1-x5);
00199     
00200     jacobian[4][0] = -timeStep*df5_dx1;
00201     jacobian[4][4] = 1-timeStep*df5_dx5;
00202     jacobian[4][5] = -timeStep*df5_dx6;
00203     
00204     // f6
00205     double df6_dx6 = mMu - 2*mMu*x6/mMstar;
00206     
00207     jacobian[5][5] = 1-timeStep*df6_dx6;
00208 }
00209 
00210 bool TysonNovak2001OdeSystem::CalculateStoppingEvent(double time, const std::vector<double> &rY)
00211 {
00212     std::vector<double> dy(rY.size());
00213     EvaluateYDerivatives(time, rY, dy);
00214     
00215     // Only call this a stopping condition if the mass of the cell is over 0.6
00216     // (normally cycles from 0.5-1.0 ish!)
00217     return ( (rY[5] > 0.6 )&& (rY[0] < mCycB_threshold) && dy[0] < 0.0 );
00218 }
00219 
00220 double TysonNovak2001OdeSystem::CalculateRootFunction(double time, const std::vector<double> &rY)
00221 {
00222     std::vector<double> dy(rY.size());
00223     EvaluateYDerivatives(time, rY, dy);
00224     
00225     // Only call this a stopping condition if the mass of the cell is over 0.6
00226     // (normally cycles from 0.5-1.0 ish!)
00227     if (rY[5]<0.6)
00228     {
00229         return 1.0;
00230     }
00231     
00232     if (dy[0] >= 0.0)
00233     {
00234         return 1.0;
00235     }
00236     return rY[0]-mCycB_threshold;
00237 }
00238 
00239 template<>
00240 void OdeSystemInformation<TysonNovak2001OdeSystem>::Initialise(void)
00241 {
00242     /*
00243      * Initialise state variables.
00244      * 
00245      * These initial conditions are the approximate steady state
00246      * solution values while the commented out conditions are taken
00247      * from the Tyson and Novak 2001 paper.
00248      */
00249     this->mVariableNames.push_back("CycB");
00250     this->mVariableUnits.push_back("nM");
00251 //    this->mInitialConditions.push_back(0.1);
00252     this->mInitialConditions.push_back(0.099999999999977);
00253     
00254     this->mVariableNames.push_back("Cdh1");
00255     this->mVariableUnits.push_back("nM");
00256 //    this->mInitialConditions.push_back(9.8770e-01);
00257     this->mInitialConditions.push_back(0.989026454281841);
00258     
00259     this->mVariableNames.push_back("Cdc20T");
00260     this->mVariableUnits.push_back("nM");
00261 //    this->mInitialConditions.push_back(1.5011e+00);
00262     this->mInitialConditions.push_back(1.547942029285891);
00263     
00264     this->mVariableNames.push_back("Cdc20A");
00265     this->mVariableUnits.push_back("nM");
00266 //    this->mInitialConditions.push_back(1.2924e+00);
00267     this->mInitialConditions.push_back(1.421110920135839);
00268     
00269     this->mVariableNames.push_back("IEP");
00270     this->mVariableUnits.push_back("nM");
00271 //    this->mInitialConditions.push_back(6.5405e-01);
00272     this->mInitialConditions.push_back(0.672838844290094);
00273     
00274     this->mVariableNames.push_back("mass");
00275     this->mVariableUnits.push_back("");
00276 //    this->mInitialConditions.push_back(4.7039e-01);
00277     this->mInitialConditions.push_back(0.970831277863956 / 2);
00278     
00279     this->mInitialised = true;
00280 }

Generated on Wed Mar 18 12:51:49 2009 for Chaste by  doxygen 1.5.5