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

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