NhsContractionModel.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 "NhsContractionModel.hpp"
00029 #include "OdeSystemInformation.hpp"
00030 #include "EulerIvpOdeSolver.hpp"
00031 #include "UblasCustomFunctions.hpp"
00032 
00033 #include <cmath>
00034 
00035 
00036 
00037 //
00038 // Model-scope constant parameters
00039 //
00040 const double NhsContractionModel::mKon = 100;
00041 const double NhsContractionModel::mKrefoff = 0.2;
00042 const double NhsContractionModel::mGamma = 2;
00043 const double NhsContractionModel::mCalciumTroponinMax = 0.07;
00044 const double NhsContractionModel::mAlphaR1 = 0.002;
00045 const double NhsContractionModel::mAlphaR2 = 0.0017;
00046 const double NhsContractionModel::mKZ = 0.15;
00047 const unsigned NhsContractionModel::mNr = 3u;
00048 const double NhsContractionModel::mBeta1 = -4;
00049 const double NhsContractionModel::mAlpha0 = 0.008;
00050 const unsigned NhsContractionModel::mN = 3u;
00051 const double NhsContractionModel::mZp = 0.85;
00052 const double NhsContractionModel::mCalcium50ref = 0.00105;
00053 const double NhsContractionModel::mTref = 56.2;
00054 const double NhsContractionModel::mBeta0 = 4.9;
00055 const double NhsContractionModel::mA = 0.35;
00056 const double NhsContractionModel::mA1 = -29;
00057 const double NhsContractionModel::mA2 = 138;
00058 const double NhsContractionModel::mA3 = 129;
00059 const double NhsContractionModel::mAlpha1 = 0.03;
00060 const double NhsContractionModel::mAlpha2 = 0.130;
00061 const double NhsContractionModel::mAlpha3 = 0.625;
00062 
00063 
00064 /*
00065  * ============================== PRIVATE FUNCTIONS =====================================
00066  */
00067 void NhsContractionModel::CalculateCalciumTrop50()
00068 {
00069     double Ca50ref_times_one_plus_beta1_times_lam_minus_one = mCalcium50ref * (1 + mBeta1*(mLambda-1));
00070     double one_plus_beta0_times_lam_minus_one_over_two_gamma = (1 + mBeta0*(mLambda-1))/(2*mGamma);
00071 
00072     mCalciumTrop50 = mCalciumTroponinMax * Ca50ref_times_one_plus_beta1_times_lam_minus_one;
00073     mCalciumTrop50 /= (Ca50ref_times_one_plus_beta1_times_lam_minus_one + (1-one_plus_beta0_times_lam_minus_one_over_two_gamma)*mKrefoff/mKon);
00074 }
00075 
00076 
00077 double NhsContractionModel::CalculateT0(double z)
00078 {
00079     double calcium_ratio_to_n = SmallPow(mCalciumTrop50/mCalciumTroponinMax, mN);
00080 
00081     double z_max = mAlpha0 - mK2*calcium_ratio_to_n;
00082     z_max /= mAlpha0 + (mAlphaR1 + mK1)*calcium_ratio_to_n;
00083 
00084     return z * mTref * (1+mBeta0*(mLambda-1)) / z_max;
00085 }
00086 
00087 
00088 
00089 /*
00090  * ============================== PUBLIC FUNCTIONS =====================================
00091  */
00092 
00093 NhsContractionModel::NhsContractionModel()
00094     :   AbstractOdeBasedContractionModel(5) // five state variables
00095 {
00096     mpSystemInfo = OdeSystemInformation<NhsContractionModel>::Instance();
00097     SetStateVariables(GetInitialConditions());
00098 
00099     mLambda = 1.0;
00100     mDLambdaDt = 0.0;
00101     mCalciumI = 0.0;
00102 
00103     // Initialise mCalciumTrop50!!
00104     CalculateCalciumTrop50();
00105 
00106     double zp_to_n_plus_K_to_n = SmallPow(mZp,mNr) + SmallPow(mKZ,mNr);
00107 
00108     mK1 = mAlphaR2 * SmallPow(mZp,mNr-1) * mNr * SmallPow(mKZ,mNr);
00109     mK1 /= zp_to_n_plus_K_to_n * zp_to_n_plus_K_to_n;
00110 
00111     mK2 = mAlphaR2 * SmallPow(mZp,mNr)/zp_to_n_plus_K_to_n;
00112     mK2 *= 1 - mNr*SmallPow(mKZ,mNr)/zp_to_n_plus_K_to_n;
00113 }
00114 
00115 void NhsContractionModel::SetStretchAndStretchRate(double lambda, double dlambdaDt)
00116 {
00117     assert(lambda>0.0);
00118     mLambda = lambda;
00119     mDLambdaDt = dlambdaDt;
00120     // lambda changed so update mCalciumTrop50!!
00121     CalculateCalciumTrop50();
00122 }
00123 
00124 void NhsContractionModel::SetInputParameters(ContractionModelInputParameters& rInputParameters)
00125 {
00126     assert(rInputParameters.intracellularCalciumConcentration != DOUBLE_UNSET);
00127     assert(rInputParameters.intracellularCalciumConcentration > 0.0);
00128     mCalciumI = rInputParameters.intracellularCalciumConcentration;
00129 }
00130 
00131 void NhsContractionModel::SetIntracellularCalciumConcentration(double calciumConcentration)
00132 {
00133     assert(calciumConcentration > 0.0);
00134     mCalciumI = calciumConcentration;
00135 }
00136 
00137 double NhsContractionModel::GetCalciumTroponinValue()
00138 {
00139     return mStateVariables[0];
00140 }
00141 
00142 void NhsContractionModel::EvaluateYDerivatives(double time,
00143                                                const std::vector<double> &rY,
00144                                                std::vector<double> &rDY)
00145 {
00147 
00148     const double& calcium_troponin = rY[0];
00149     const double& z = rY[1];
00150     const double& Q1 = rY[2];
00151     const double& Q2 = rY[3];
00152     const double& Q3 = rY[4];
00153 
00154     // check the state vars are in the expected range
00155     #define COVERAGE_IGNORE
00156     if(calcium_troponin < 0)
00157     {
00158         EXCEPTION("CalciumTrop concentration went negative");
00159     }
00160     if(z<0)
00161     {
00162         EXCEPTION("z went negative");
00163     }
00164     if(z>1)
00165     {
00166         EXCEPTION("z became greater than 1");
00167     }
00168     #undef COVERAGE_IGNORE
00169 
00170 
00171     double Q = Q1 + Q2 + Q3;
00172     double T0 = CalculateT0(z);
00173 
00174     double Ta;
00175     if(Q>0)
00176     {
00177         Ta = T0*(1+(2+mA)*Q)/(1+Q);
00178     }
00179     else
00180     {
00181         Ta = T0*(1+mA*Q)/(1-Q);
00182     }
00183 
00184     rDY[0] =   mKon * mCalciumI * ( mCalciumTroponinMax - calcium_troponin)
00185              - mKrefoff * (1-Ta/(mGamma*mTref)) * calcium_troponin;
00186 
00187     double ca_trop_to_ca_trop50_ratio_to_n = SmallPow(calcium_troponin/mCalciumTrop50, mN);
00188 
00189     rDY[1] =   mAlpha0 * ca_trop_to_ca_trop50_ratio_to_n * (1-z)
00190              - mAlphaR1 * z
00191              - mAlphaR2 * SmallPow(z,mNr) / (SmallPow(z,mNr) + SmallPow(mKZ,mNr));
00192 
00193 
00194     rDY[2] = mA1 * mDLambdaDt - mAlpha1 * Q1;
00195     rDY[3] = mA2 * mDLambdaDt - mAlpha2 * Q2;
00196     rDY[4] = mA3 * mDLambdaDt - mAlpha3 * Q3;
00197 }
00198 
00199 
00200 double NhsContractionModel::GetActiveTension()
00201 {
00202     double T0 = CalculateT0(mStateVariables[1]);
00203     double Q = mStateVariables[2]+mStateVariables[3]+mStateVariables[4];
00204 
00205     if(Q>0)
00206     {
00207         return T0*(1+(2+mA)*Q)/(1+Q);
00208     }
00209     else
00210     {
00211         return T0*(1+mA*Q)/(1-Q);
00212     }
00213 }
00214 
00215 template<>
00216 void OdeSystemInformation<NhsContractionModel>::Initialise(void)
00217 {
00218     this->mVariableNames.push_back("CalciumTroponin");
00219     this->mVariableUnits.push_back("microMols");
00220     this->mInitialConditions.push_back(0);
00221 
00222     this->mVariableNames.push_back("z");
00223     this->mVariableUnits.push_back("");
00224     this->mInitialConditions.push_back(0);
00225 
00226     this->mVariableNames.push_back("Q1");
00227     this->mVariableUnits.push_back("");
00228     this->mInitialConditions.push_back(0);
00229 
00230     this->mVariableNames.push_back("Q2");
00231     this->mVariableUnits.push_back("");
00232     this->mInitialConditions.push_back(0);
00233 
00234     this->mVariableNames.push_back("Q3");
00235     this->mVariableUnits.push_back("");
00236     this->mInitialConditions.push_back(0);
00237 
00238     this->mInitialised = true;
00239 }

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