LuoRudyIModel1991OdeSystem.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2010
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 "LuoRudyIModel1991OdeSystem.hpp"
00029 #include "OdeSystemInformation.hpp"
00030 #include <cmath>
00031 //#include <iostream>
00032 #include "Exception.hpp"
00033 
00034 
00035 //
00036 // Model-scope constant parameters
00037 //
00038 const double LuoRudyIModel1991OdeSystem::membrane_C = 1.0;
00039 const double LuoRudyIModel1991OdeSystem::membrane_F = 96484.6;
00040 const double LuoRudyIModel1991OdeSystem::membrane_R = 8314;
00041 const double LuoRudyIModel1991OdeSystem::membrane_T = 310.0;
00042 const double LuoRudyIModel1991OdeSystem::background_current_E_b = -59.87;
00043 const double LuoRudyIModel1991OdeSystem::background_current_g_b = 0.03921;
00044 const double LuoRudyIModel1991OdeSystem::fast_sodium_current_g_Na = 23.0;
00045 const double LuoRudyIModel1991OdeSystem::ionic_concentrations_Ki = 145.0;
00046 const double LuoRudyIModel1991OdeSystem::ionic_concentrations_Ko = 5.4;
00047 const double LuoRudyIModel1991OdeSystem::ionic_concentrations_Nai = 18.0;
00048 const double LuoRudyIModel1991OdeSystem::ionic_concentrations_Nao = 140.0;
00049 const double LuoRudyIModel1991OdeSystem::plateau_potassium_current_g_Kp = 0.0183;
00050 const double LuoRudyIModel1991OdeSystem::time_dependent_potassium_current_PR_NaK = 0.01833;
00051 
00052 
00053 
00057 LuoRudyIModel1991OdeSystem::LuoRudyIModel1991OdeSystem(
00058     boost::shared_ptr<AbstractIvpOdeSolver> pSolver,
00059     boost::shared_ptr<AbstractStimulusFunction> pIntracellularStimulus)
00060         : AbstractCardiacCell(pSolver, 8, 4, pIntracellularStimulus)
00061 {
00062     mpSystemInfo = OdeSystemInformation<LuoRudyIModel1991OdeSystem>::Instance();
00063 
00064     // set the final paramter
00065     fast_sodium_current_E_Na = ((membrane_R * membrane_T) / membrane_F) *
00066                                log(ionic_concentrations_Nao / ionic_concentrations_Nai);
00067 
00068     Init();
00069 }
00070 
00074 LuoRudyIModel1991OdeSystem::~LuoRudyIModel1991OdeSystem(void)
00075 {
00076 }
00077 
00078 double LuoRudyIModel1991OdeSystem::GetIntracellularCalciumConcentration()
00079 {
00080     return mStateVariables[3];
00081 }
00082 
00092 void LuoRudyIModel1991OdeSystem::EvaluateYDerivatives(double time,
00093                                                       const std::vector<double> &rY,
00094                                                       std::vector<double> &rDY)
00095 {
00096     double fast_sodium_current_h_gate_h = rY[0];
00097     double fast_sodium_current_j_gate_j = rY[1];
00098     double fast_sodium_current_m_gate_m = rY[2];
00099     double intracellular_calcium_concentration_Cai = rY[3];
00100     double membrane_V = rY[4];
00101     double slow_inward_current_d_gate_d = rY[5];
00102     double slow_inward_current_f_gate_f = rY[6];
00103     double time_dependent_potassium_current_X_gate_X = rY[7];
00104 
00105     double background_current_i_b = background_current_g_b*(membrane_V-background_current_E_b);
00106 
00107     double fast_sodium_current_h_gate_alpha_h;
00108 
00109     if (membrane_V < -40.0)
00110     {
00111         fast_sodium_current_h_gate_alpha_h = 0.135*exp((80.0+membrane_V)/-6.8);
00112     }
00113     else
00114     {
00115         fast_sodium_current_h_gate_alpha_h = 0.0;
00116     }
00117 
00118     double fast_sodium_current_h_gate_beta_h;
00119 
00120     if (membrane_V < -40.0)
00121     {
00122         fast_sodium_current_h_gate_beta_h = 3.56*exp(0.079*membrane_V)+3.1e5*exp(0.35*membrane_V);
00123     }
00124     else
00125     {
00126         fast_sodium_current_h_gate_beta_h = 1.0/(0.13*(1.0+exp((membrane_V+10.66)/-11.1)));
00127     }
00128 
00129     double fast_sodium_current_h_gate_h_prime = fast_sodium_current_h_gate_alpha_h*(1.0-fast_sodium_current_h_gate_h)-fast_sodium_current_h_gate_beta_h*fast_sodium_current_h_gate_h;
00130 
00131     double fast_sodium_current_j_gate_alpha_j;
00132 
00133     if (membrane_V < -40.0)
00134     {
00135         fast_sodium_current_j_gate_alpha_j = (-1.2714e5*exp(0.2444*membrane_V)-3.474e-5*exp(-0.04391*membrane_V))*(membrane_V+37.78)/(1.0+exp(0.311*(membrane_V+79.23)));
00136     }
00137     else
00138     {
00139         fast_sodium_current_j_gate_alpha_j = 0.0;
00140     }
00141 
00142     double fast_sodium_current_j_gate_beta_j;
00143 
00144     if (membrane_V < -40.0)
00145     {
00146         fast_sodium_current_j_gate_beta_j = 0.1212*exp(-0.01052*membrane_V)/(1.0+exp(-0.1378*(membrane_V+40.14)));
00147     }
00148     else
00149     {
00150         fast_sodium_current_j_gate_beta_j = 0.3*exp(-2.535e-7*membrane_V)/(1.0+exp(-0.1*(membrane_V+32.0)));
00151     }
00152 
00153     double fast_sodium_current_j_gate_j_prime = fast_sodium_current_j_gate_alpha_j*(1.0-fast_sodium_current_j_gate_j)-fast_sodium_current_j_gate_beta_j*fast_sodium_current_j_gate_j;
00154     double fast_sodium_current_m_gate_alpha_m = 0.32*(membrane_V+47.13)/(1.0-exp(-0.1*(membrane_V+47.13)));
00155     double fast_sodium_current_m_gate_beta_m = 0.08*exp(-membrane_V/11.0);
00156     double fast_sodium_current_m_gate_m_prime = fast_sodium_current_m_gate_alpha_m*(1.0-fast_sodium_current_m_gate_m)-fast_sodium_current_m_gate_beta_m*fast_sodium_current_m_gate_m;
00157     double fast_sodium_current_i_Na = fast_sodium_current_g_Na*pow(fast_sodium_current_m_gate_m, 3.0)*fast_sodium_current_h_gate_h*fast_sodium_current_j_gate_j*(membrane_V-fast_sodium_current_E_Na);
00158 
00159     double slow_inward_current_d_gate_alpha_d = 0.095*exp(-0.01*(membrane_V-5.0))/(1.0+exp(-0.072*(membrane_V-5.0)));
00160     double slow_inward_current_d_gate_beta_d = 0.07*exp(-0.017*(membrane_V+44.0))/(1.0+exp(0.05*(membrane_V+44.0)));
00161     double slow_inward_current_d_gate_d_prime = slow_inward_current_d_gate_alpha_d*(1.0-slow_inward_current_d_gate_d)-slow_inward_current_d_gate_beta_d*slow_inward_current_d_gate_d;
00162 
00163     double slow_inward_current_f_gate_alpha_f = 0.012*exp(-0.008*(membrane_V+28.0))/(1.0+exp(0.15*(membrane_V+28.0)));
00164     double slow_inward_current_f_gate_beta_f = 0.0065*exp(-0.02*(membrane_V+30.0))/(1.0+exp(-0.2*(membrane_V+30.0)));
00165     double slow_inward_current_f_gate_f_prime = slow_inward_current_f_gate_alpha_f*(1.0-slow_inward_current_f_gate_f)-slow_inward_current_f_gate_beta_f*slow_inward_current_f_gate_f;
00166 
00167     double slow_inward_current_E_si = 7.7-13.0287*log(intracellular_calcium_concentration_Cai);
00168     double slow_inward_current_i_si = 0.09*slow_inward_current_d_gate_d*slow_inward_current_f_gate_f*(membrane_V-slow_inward_current_E_si);
00169     double intracellular_calcium_concentration_Cai_prime = -1e-4*slow_inward_current_i_si+0.07*(1e-4-intracellular_calcium_concentration_Cai);
00170     double time_dependent_potassium_current_g_K = 0.282*sqrt(ionic_concentrations_Ko/5.4);
00171 
00172     double time_dependent_potassium_current_Xi_gate_Xi;
00173 
00174     if (membrane_V > -100.0)
00175     {
00176         time_dependent_potassium_current_Xi_gate_Xi = 2.837*(exp(0.04*(membrane_V+77.0))-1.0)/((membrane_V+77.0)*exp(0.04*(membrane_V+35.0)));
00177     }
00178     else
00179     {
00180         #define COVERAGE_IGNORE
00181         time_dependent_potassium_current_Xi_gate_Xi = 1.0;
00182         #undef COVERAGE_IGNORE
00183     }
00184 
00185     double time_dependent_potassium_current_X_gate_alpha_X = 0.0005*exp(0.083*(membrane_V+50.0))/(1.0+exp(0.057*(membrane_V+50.0)));
00186     double time_dependent_potassium_current_X_gate_beta_X = 0.0013*exp(-0.06*(membrane_V+20.0))/(1.0+exp(-0.04*(membrane_V+20.0)));
00187     double time_dependent_potassium_current_X_gate_X_prime = time_dependent_potassium_current_X_gate_alpha_X*(1.0-time_dependent_potassium_current_X_gate_X)-time_dependent_potassium_current_X_gate_beta_X*time_dependent_potassium_current_X_gate_X;
00188 
00189     double time_dependent_potassium_current_E_K = ((membrane_R*membrane_T)/membrane_F)*log((ionic_concentrations_Ko+time_dependent_potassium_current_PR_NaK*ionic_concentrations_Nao)/(ionic_concentrations_Ki+time_dependent_potassium_current_PR_NaK*ionic_concentrations_Nai));
00190     double time_dependent_potassium_current_i_K = time_dependent_potassium_current_g_K*time_dependent_potassium_current_X_gate_X*time_dependent_potassium_current_Xi_gate_Xi*(membrane_V-time_dependent_potassium_current_E_K);
00191     double time_independent_potassium_current_g_K1 = 0.6047*sqrt(ionic_concentrations_Ko/5.4);
00192     double time_independent_potassium_current_E_K1 =((membrane_R*membrane_T)/membrane_F)*log(ionic_concentrations_Ko/ionic_concentrations_Ki);
00193     double time_independent_potassium_current_K1_gate_alpha_K1 = 1.02/(1.0+exp(0.2385*(membrane_V-time_independent_potassium_current_E_K1-59.215)));
00194     double time_independent_potassium_current_K1_gate_beta_K1 = (0.49124*exp(0.08032*(membrane_V+5.476-time_independent_potassium_current_E_K1))+exp(0.06175*(membrane_V-(time_independent_potassium_current_E_K1+594.31))))/(1.0+exp(-0.5143*(membrane_V-time_independent_potassium_current_E_K1+4.753)));
00195     double time_independent_potassium_current_K1_gate_K1_infinity = time_independent_potassium_current_K1_gate_alpha_K1/(time_independent_potassium_current_K1_gate_alpha_K1+time_independent_potassium_current_K1_gate_beta_K1);
00196     double time_independent_potassium_current_i_K1 = time_independent_potassium_current_g_K1*time_independent_potassium_current_K1_gate_K1_infinity*(membrane_V-time_independent_potassium_current_E_K1);
00197     double plateau_potassium_current_Kp = 1.0/(1.0+exp((7.488-membrane_V)/5.98));
00198     double plateau_potassium_current_E_Kp = time_independent_potassium_current_E_K1;
00199     double plateau_potassium_current_i_Kp = plateau_potassium_current_g_Kp*plateau_potassium_current_Kp*(membrane_V-plateau_potassium_current_E_Kp);
00200     double i_stim = GetStimulus(time);
00201 
00202     //calculate dV
00203     double membrane_V_prime = (-1.0/membrane_C)*(fast_sodium_current_i_Na+slow_inward_current_i_si+time_dependent_potassium_current_i_K+time_independent_potassium_current_i_K1+plateau_potassium_current_i_Kp+background_current_i_b + i_stim);
00204     // do not update voltage if the mSetVoltageDerivativeToZero flag has been set
00205     if (mSetVoltageDerivativeToZero)
00206     {
00207         membrane_V_prime = 0;
00208     }
00209 
00210     rDY[0] = fast_sodium_current_h_gate_h_prime;
00211     rDY[1] = fast_sodium_current_j_gate_j_prime;
00212     rDY[2] = fast_sodium_current_m_gate_m_prime;
00213     rDY[3] = intracellular_calcium_concentration_Cai_prime;
00214     rDY[4] = membrane_V_prime;
00215     rDY[5] = slow_inward_current_d_gate_d_prime;
00216     rDY[6] = slow_inward_current_f_gate_f_prime;
00217     rDY[7] = time_dependent_potassium_current_X_gate_X_prime;
00218 }
00219 
00220 
00221 double LuoRudyIModel1991OdeSystem::GetIIonic()
00222 {
00223     double fast_sodium_current_h_gate_h = mStateVariables[0];
00224     double fast_sodium_current_j_gate_j = mStateVariables[1];
00225     double fast_sodium_current_m_gate_m = mStateVariables[2];
00226     double intracellular_calcium_concentration_Cai = mStateVariables[3];
00227     double membrane_V = mStateVariables[4];
00228     double slow_inward_current_d_gate_d = mStateVariables[5];
00229     double slow_inward_current_f_gate_f = mStateVariables[6];
00230     double time_dependent_potassium_current_X_gate_X = mStateVariables[7];
00231 
00232     /*
00233      * Compute the LuoRudyIModel1991OdeSystem model
00234      */
00235     double background_current_i_b = background_current_g_b*(membrane_V-background_current_E_b);
00236 
00237     double fast_sodium_current_i_Na = fast_sodium_current_g_Na*pow(fast_sodium_current_m_gate_m, 3.0)*fast_sodium_current_h_gate_h*fast_sodium_current_j_gate_j*(membrane_V-fast_sodium_current_E_Na);
00238 
00239     double slow_inward_current_E_si = 7.7-13.0287*log(intracellular_calcium_concentration_Cai);
00240     double slow_inward_current_i_si = 0.09*slow_inward_current_d_gate_d*slow_inward_current_f_gate_f*(membrane_V-slow_inward_current_E_si);
00241 
00242     double time_dependent_potassium_current_g_K = 0.282*sqrt(ionic_concentrations_Ko/5.4);
00243     double time_dependent_potassium_current_Xi_gate_Xi;
00244 
00245     // Although the equation below looks strange (particularly the arguments of the
00246     // exponentials, it is in fact correct.
00247     if (membrane_V > -100.0)
00248     {
00249         time_dependent_potassium_current_Xi_gate_Xi = 2.837*(exp(0.04*(membrane_V+77.0))-1.0)/((membrane_V+77.0)*exp(0.04*(membrane_V+35.0)));
00250     }
00251     else
00252     {
00253         #define COVERAGE_IGNORE
00254         time_dependent_potassium_current_Xi_gate_Xi = 1.0;
00255         #undef COVERAGE_IGNORE
00256     }
00257     double time_dependent_potassium_current_E_K = ((membrane_R*membrane_T)/membrane_F)*log((ionic_concentrations_Ko+time_dependent_potassium_current_PR_NaK*ionic_concentrations_Nao)/(ionic_concentrations_Ki+time_dependent_potassium_current_PR_NaK*ionic_concentrations_Nai));
00258     double time_dependent_potassium_current_i_K = time_dependent_potassium_current_g_K*time_dependent_potassium_current_X_gate_X*time_dependent_potassium_current_Xi_gate_Xi*(membrane_V-time_dependent_potassium_current_E_K);
00259 
00260     double time_independent_potassium_current_g_K1 = 0.6047*sqrt(ionic_concentrations_Ko/5.4);
00261     double time_independent_potassium_current_E_K1 =((membrane_R*membrane_T)/membrane_F)*log(ionic_concentrations_Ko/ionic_concentrations_Ki);
00262     double time_independent_potassium_current_K1_gate_alpha_K1 = 1.02/(1.0+exp(0.2385*(membrane_V-time_independent_potassium_current_E_K1-59.215)));
00263     double time_independent_potassium_current_K1_gate_beta_K1 = (0.49124*exp(0.08032*(membrane_V+5.476-time_independent_potassium_current_E_K1))+exp(0.06175*(membrane_V-(time_independent_potassium_current_E_K1+594.31))))/(1.0+exp(-0.5143*(membrane_V-time_independent_potassium_current_E_K1+4.753)));
00264     double time_independent_potassium_current_K1_gate_K1_infinity = time_independent_potassium_current_K1_gate_alpha_K1/(time_independent_potassium_current_K1_gate_alpha_K1+time_independent_potassium_current_K1_gate_beta_K1);
00265     double time_independent_potassium_current_i_K1 = time_independent_potassium_current_g_K1*time_independent_potassium_current_K1_gate_K1_infinity*(membrane_V-time_independent_potassium_current_E_K1);
00266 
00267     double plateau_potassium_current_Kp = 1.0/(1.0+exp((7.488-membrane_V)/5.98));
00268     double plateau_potassium_current_E_Kp = time_independent_potassium_current_E_K1;
00269     double plateau_potassium_current_i_Kp = plateau_potassium_current_g_Kp*plateau_potassium_current_Kp*(membrane_V-plateau_potassium_current_E_Kp);
00270 
00271     double i_ionic = fast_sodium_current_i_Na+slow_inward_current_i_si+time_dependent_potassium_current_i_K+time_independent_potassium_current_i_K1+plateau_potassium_current_i_Kp+background_current_i_b;
00272 
00273     assert(!std::isnan(i_ionic));
00274     return i_ionic;
00275 }
00276 
00277 void LuoRudyIModel1991OdeSystem::VerifyStateVariables()
00278 {
00279     const std::vector<double>& rY = rGetStateVariables();
00280 
00281     const double fast_sodium_current_h_gate_h = rY[0];            // gating
00282     const double fast_sodium_current_j_gate_j = rY[1];            // gating
00283     const double fast_sodium_current_m_gate_m = rY[2];            // gating
00284     const double intracellular_calcium_concentration_Cai = rY[3]; // concentration
00285     const double slow_inward_current_d_gate_d = rY[5];            // gating
00286     const double slow_inward_current_f_gate_f = rY[6];            // gating
00287     const double time_dependent_potassium_current_X_gate_X = rY[7]; // gating
00288 
00289     #define COVERAGE_IGNORE
00290     if (!(0.0<=fast_sodium_current_h_gate_h && fast_sodium_current_h_gate_h<=1.0))
00291     {
00292         EXCEPTION(DumpState("h gate for fast sodium current has gone out of range. Check model parameters, for example spatial stepsize"));
00293     }
00294 
00295     if (!(0.0<=fast_sodium_current_j_gate_j && fast_sodium_current_j_gate_j<=1.0))
00296     {
00297         EXCEPTION(DumpState("j gate for fast sodium current has gone out of range. Check model parameters, for example spatial stepsize"));
00298     }
00299 
00300     if (!(0.0<=fast_sodium_current_m_gate_m && fast_sodium_current_m_gate_m<=1.0))
00301     {
00302         EXCEPTION(DumpState("m gate for fast sodium current has gone out of range. Check model parameters, for example spatial stepsize"));
00303     }
00304 
00305     if (!(0.0<intracellular_calcium_concentration_Cai))
00306     {
00307         EXCEPTION(DumpState("intracellular_calcium_concentration_Cai has become non-positive, ie gone out of range. Check model parameters, for example spatial stepsize"));
00308     }
00309 
00310     if (!(0.0<=slow_inward_current_d_gate_d && slow_inward_current_d_gate_d<=1.0))
00311     {
00312         EXCEPTION(DumpState("d gate for slow inward current has gone out of range. Check model parameters, for example spatial stepsize"));
00313     }
00314 
00315     if (!(0.0<=slow_inward_current_f_gate_f && slow_inward_current_f_gate_f<=1.0))
00316     {
00317         EXCEPTION(DumpState("f gate for slow inward current has gone out of range. Check model parameters, for example spatial stepsize"));
00318     }
00319 
00320     if (!(0.0<=time_dependent_potassium_current_X_gate_X && time_dependent_potassium_current_X_gate_X<=1.0))
00321     {
00322         EXCEPTION(DumpState("X gate for time dependent potassium current has gone out of range. Check model parameters, for example spatial stepsize"));
00323     }
00324     #undef COVERAGE_IGNORE
00325 }
00326 
00327 
00328 
00329 
00330 template<>
00331 void OdeSystemInformation<LuoRudyIModel1991OdeSystem>::Initialise(void)
00332 {
00333     // State variables
00334     this->mVariableNames.push_back("h");
00335     this->mVariableUnits.push_back("");
00336     this->mInitialConditions.push_back(0.9804713);
00337 
00338     this->mVariableNames.push_back("j");
00339     this->mVariableUnits.push_back("");
00340     this->mInitialConditions.push_back(0.98767124);
00341 
00342     this->mVariableNames.push_back("m");
00343     this->mVariableUnits.push_back("");
00344     this->mInitialConditions.push_back(0.00187018);
00345 
00346     this->mVariableNames.push_back("CaI");
00347     this->mVariableUnits.push_back("mMol");
00348     this->mInitialConditions.push_back(0.0002);
00349 
00350     this->mVariableNames.push_back("V");
00351     this->mVariableUnits.push_back("mV");
00352     this->mInitialConditions.push_back(-83.853);
00353 
00354     this->mVariableNames.push_back("d");
00355     this->mVariableUnits.push_back("");
00356     this->mInitialConditions.push_back(0.00316354);
00357 
00358     this->mVariableNames.push_back("f");
00359     this->mVariableUnits.push_back("");
00360     this->mInitialConditions.push_back(0.99427859);
00361 
00362     this->mVariableNames.push_back("x");
00363     this->mVariableUnits.push_back("");
00364     this->mInitialConditions.push_back(0.16647703);
00365 
00366     this->mInitialised = true;
00367 }
00368 
00369 
00370 // Serialization for Boost >= 1.36
00371 #include "SerializationExportWrapperForCpp.hpp"
00372 CHASTE_CLASS_EXPORT(LuoRudyIModel1991OdeSystem);

Generated by  doxygen 1.6.2