BackwardEulerLuoRudyIModel1991.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 "BackwardEulerLuoRudyIModel1991.hpp"
00029 #include "Exception.hpp"
00030 #include "CardiacNewtonSolver.hpp"
00031 #include "OdeSystemInformation.hpp"
00032 
00033 #include <cmath>
00034 #include <cassert>
00035 //#include <iostream>
00036 
00037 
00038 //
00039 // Constant model-scope parameters
00040 //
00041 const double BackwardEulerLuoRudyIModel1991::membrane_C = 1.0;
00042 const double BackwardEulerLuoRudyIModel1991::membrane_F = 96484.6;
00043 const double BackwardEulerLuoRudyIModel1991::membrane_R = 8314;
00044 const double BackwardEulerLuoRudyIModel1991::membrane_T = 310.0;
00045 const double BackwardEulerLuoRudyIModel1991::background_current_E_b = -59.87;
00046 const double BackwardEulerLuoRudyIModel1991::background_current_g_b = 0.03921;
00047 const double BackwardEulerLuoRudyIModel1991::fast_sodium_current_g_Na = 23.0;
00048 const double BackwardEulerLuoRudyIModel1991::ionic_concentrations_Ki = 145.0;
00049 const double BackwardEulerLuoRudyIModel1991::ionic_concentrations_Ko = 5.4;
00050 const double BackwardEulerLuoRudyIModel1991::ionic_concentrations_Nai = 18.0;
00051 const double BackwardEulerLuoRudyIModel1991::ionic_concentrations_Nao = 140.0;
00052 const double BackwardEulerLuoRudyIModel1991::plateau_potassium_current_g_Kp = 0.0183;
00053 const double BackwardEulerLuoRudyIModel1991::time_dependent_potassium_current_PR_NaK = 0.01833;
00054 
00055 
00056 BackwardEulerLuoRudyIModel1991::BackwardEulerLuoRudyIModel1991(
00057     boost::shared_ptr<AbstractStimulusFunction> pIntracellularStimulus)
00058         : AbstractBackwardEulerCardiacCell<1>(8, 4, pIntracellularStimulus)
00059 {
00060     Init();
00061 }
00062 
00063 BackwardEulerLuoRudyIModel1991::BackwardEulerLuoRudyIModel1991(
00064     boost::shared_ptr<AbstractIvpOdeSolver> /* unused */,
00065     boost::shared_ptr<AbstractStimulusFunction> pIntracellularStimulus)
00066         : AbstractBackwardEulerCardiacCell<1>(8, 4, pIntracellularStimulus)
00067 {
00068     Init();
00069 }
00070 
00071 BackwardEulerLuoRudyIModel1991::~BackwardEulerLuoRudyIModel1991(void)
00072 {
00073 }
00074 
00075 
00076 void BackwardEulerLuoRudyIModel1991::Init()
00077 {
00078     mpSystemInfo = OdeSystemInformation<BackwardEulerLuoRudyIModel1991>::Instance();
00079 
00080     // set final parameter member variable
00081     fast_sodium_current_E_Na = ((membrane_R * membrane_T) / membrane_F) *
00082                                 log(ionic_concentrations_Nao / ionic_concentrations_Nai);
00083 
00084     AbstractCardiacCell::Init();
00085 }
00086 
00087 template<>
00088 void OdeSystemInformation<BackwardEulerLuoRudyIModel1991>::Initialise(void)
00089 {
00090     this->mVariableNames.push_back("h");
00091     this->mVariableUnits.push_back("");
00092     this->mInitialConditions.push_back(0.9804713);
00093 
00094     this->mVariableNames.push_back("j");
00095     this->mVariableUnits.push_back("");
00096     this->mInitialConditions.push_back(0.98767124);
00097 
00098     this->mVariableNames.push_back("m");
00099     this->mVariableUnits.push_back("");
00100     this->mInitialConditions.push_back(0.00187018);
00101 
00102     this->mVariableNames.push_back("CaI");
00103     this->mVariableUnits.push_back("mMol");
00104     this->mInitialConditions.push_back(0.0002);
00105 
00106     this->mVariableNames.push_back("V");
00107     this->mVariableUnits.push_back("mV");
00108     this->mInitialConditions.push_back(-83.853);
00109 
00110     this->mVariableNames.push_back("d");
00111     this->mVariableUnits.push_back("");
00112     this->mInitialConditions.push_back(0.00316354);
00113 
00114     this->mVariableNames.push_back("f");
00115     this->mVariableUnits.push_back("");
00116     this->mInitialConditions.push_back(0.99427859);
00117 
00118     this->mVariableNames.push_back("x");
00119     this->mVariableUnits.push_back("");
00120     this->mInitialConditions.push_back(0.16647703);
00121 
00122     this->mInitialised = true;
00123 }
00124 
00125 double BackwardEulerLuoRudyIModel1991::GetIntracellularCalciumConcentration()
00126 {
00127     return mStateVariables[3];
00128 }
00129 
00130 
00131 void BackwardEulerLuoRudyIModel1991::UpdateTransmembranePotential(double time)
00132 {
00133     // Compute next value of V
00134     std::vector<double> &rY = rGetStateVariables();
00135 
00136     double fast_sodium_current_h_gate_h = rY[0];
00137     double fast_sodium_current_j_gate_j = rY[1];
00138     double fast_sodium_current_m_gate_m = rY[2];
00139     double intracellular_calcium_concentration_Cai = rY[3];
00140     double membrane_V = rY[4];
00141     double slow_inward_current_d_gate_d = rY[5];
00142     double slow_inward_current_f_gate_f = rY[6];
00143     double time_dependent_potassium_current_X_gate_X = rY[7];
00144 
00145     double background_current_i_b = background_current_g_b*(membrane_V-background_current_E_b);
00146     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);
00147     double slow_inward_current_E_si = 7.7-13.0287*log(intracellular_calcium_concentration_Cai);
00148     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);
00149 
00150     double time_dependent_potassium_current_g_K = 0.282*sqrt(ionic_concentrations_Ko/5.4);
00151     double time_dependent_potassium_current_Xi_gate_Xi;
00152 
00153     if (membrane_V > -100.0)
00154     {
00155         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)));
00156     }
00157     else
00158     {
00159         #define COVERAGE_IGNORE
00160         time_dependent_potassium_current_Xi_gate_Xi = 1.0;
00161         #undef COVERAGE_IGNORE
00162     }
00163 
00164     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));
00165     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);
00166     double time_independent_potassium_current_g_K1 = 0.6047*sqrt(ionic_concentrations_Ko/5.4);
00167     double time_independent_potassium_current_E_K1 =((membrane_R*membrane_T)/membrane_F)*log(ionic_concentrations_Ko/ionic_concentrations_Ki);
00168     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)));
00169     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)));
00170     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);
00171     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);
00172     double plateau_potassium_current_Kp = 1.0/(1.0+exp((7.488-membrane_V)/5.98));
00173     double plateau_potassium_current_E_Kp = time_independent_potassium_current_E_K1;
00174     double plateau_potassium_current_i_Kp = plateau_potassium_current_g_Kp*plateau_potassium_current_Kp*(membrane_V-plateau_potassium_current_E_Kp);
00175     double i_stim = GetStimulus(time);
00176 
00177     //calculate dV
00178     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);
00179 
00180     rY[4] += mDt * membrane_V_prime;
00181 }
00182 
00183 void BackwardEulerLuoRudyIModel1991::ComputeOneStepExceptVoltage(double tStart)
00184 {
00185     // This method computes one timestep for all state variables except V, using
00186     // backward Euler.
00187     std::vector<double>& rY = rGetStateVariables();
00188     double fast_sodium_current_h_gate_h = rY[0];
00189     double fast_sodium_current_j_gate_j = rY[1];
00190     double fast_sodium_current_m_gate_m = rY[2];
00191 
00192     double membrane_V = rY[4];
00193     double slow_inward_current_d_gate_d = rY[5];
00194     double slow_inward_current_f_gate_f = rY[6];
00195     double time_dependent_potassium_current_X_gate_X = rY[7];
00196 
00197     double fast_sodium_current_h_gate_alpha_h;
00198 
00199     if (membrane_V < -40.0)
00200     {
00201         fast_sodium_current_h_gate_alpha_h = 0.135*exp((80.0+membrane_V)/-6.8);
00202     }
00203     else
00204     {
00205         fast_sodium_current_h_gate_alpha_h = 0.0;
00206     }
00207 
00208     double fast_sodium_current_h_gate_beta_h;
00209 
00210     if (membrane_V < -40.0)
00211     {
00212         fast_sodium_current_h_gate_beta_h = 3.56*exp(0.079*membrane_V)+3.1e5*exp(0.35*membrane_V);
00213     }
00214     else
00215     {
00216         fast_sodium_current_h_gate_beta_h = 1.0/(0.13*(1.0+exp((membrane_V+10.66)/-11.1)));
00217     }
00218 
00219     double fast_sodium_current_j_gate_alpha_j;
00220 
00221     if (membrane_V < -40.0)
00222     {
00223         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)));
00224     }
00225     else
00226     {
00227         fast_sodium_current_j_gate_alpha_j = 0.0;
00228     }
00229 
00230     double fast_sodium_current_j_gate_beta_j;
00231 
00232     if (membrane_V < -40.0)
00233     {
00234         fast_sodium_current_j_gate_beta_j = 0.1212*exp(-0.01052*membrane_V)/(1.0+exp(-0.1378*(membrane_V+40.14)));
00235     }
00236     else
00237     {
00238         fast_sodium_current_j_gate_beta_j = 0.3*exp(-2.535e-7*membrane_V)/(1.0+exp(-0.1*(membrane_V+32.0)));
00239     }
00240 
00241     double fast_sodium_current_m_gate_alpha_m = 0.32*(membrane_V+47.13)/(1.0-exp(-0.1*(membrane_V+47.13)));
00242     double fast_sodium_current_m_gate_beta_m = 0.08*exp(-membrane_V/11.0);
00243 
00244     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)));
00245     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)));
00246 
00247     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)));
00248     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)));
00249 
00250     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)));
00251     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)));
00252 
00253     // update h
00254     fast_sodium_current_h_gate_h =   (fast_sodium_current_h_gate_h + fast_sodium_current_h_gate_alpha_h*mDt)
00255                                      / (1 + (fast_sodium_current_h_gate_alpha_h + fast_sodium_current_h_gate_beta_h)*mDt);
00256 
00257     // update m
00258     fast_sodium_current_m_gate_m =   (fast_sodium_current_m_gate_m + fast_sodium_current_m_gate_alpha_m*mDt)
00259                                      / (1 + (fast_sodium_current_m_gate_alpha_m + fast_sodium_current_m_gate_beta_m)*mDt);
00260 
00261     // update j
00262     fast_sodium_current_j_gate_j =   (fast_sodium_current_j_gate_j + fast_sodium_current_j_gate_alpha_j*mDt)
00263                                      / (1 + (fast_sodium_current_j_gate_alpha_j + fast_sodium_current_j_gate_beta_j)*mDt);
00264 
00265     // update d
00266     slow_inward_current_d_gate_d =   (slow_inward_current_d_gate_d + slow_inward_current_d_gate_alpha_d*mDt)
00267                                      / (1 + (slow_inward_current_d_gate_alpha_d + slow_inward_current_d_gate_beta_d)*mDt);
00268 
00269     // update f
00270     slow_inward_current_f_gate_f =   (slow_inward_current_f_gate_f + slow_inward_current_f_gate_alpha_f*mDt)
00271                                      / (1 + (slow_inward_current_f_gate_alpha_f + slow_inward_current_f_gate_beta_f)*mDt);
00272 
00273     // update X
00274     time_dependent_potassium_current_X_gate_X =   ( time_dependent_potassium_current_X_gate_X + time_dependent_potassium_current_X_gate_alpha_X*mDt)
00275                                                   / ( 1 + (time_dependent_potassium_current_X_gate_alpha_X+time_dependent_potassium_current_X_gate_beta_X)*mDt);
00276 
00277     rY[0] =  fast_sodium_current_h_gate_h;
00278     rY[1] =  fast_sodium_current_j_gate_j;
00279     rY[2] =  fast_sodium_current_m_gate_m;
00280 
00281     rY[5] =  slow_inward_current_d_gate_d;
00282     rY[6] =  slow_inward_current_f_gate_f;
00283     rY[7] =  time_dependent_potassium_current_X_gate_X;
00284 
00285     // finally, use the newton method to calculate the Calcium concentration at the next time
00286     double guess[1] = {rY[3]};
00287 
00288     CardiacNewtonSolver<1> *solver = CardiacNewtonSolver<1>::Instance();
00289     solver->Solve(*this, tStart, guess);
00290 
00291     rY[3] = guess[0];
00292 
00293     // Timings:
00294     //  Hardcoded one-variable newton:
00295     //    Forward: 1.66
00296     //    Backward: 1.29
00297     //    Backward (long dt): 0.03
00298     //  Original DoNewton:
00299     //    Forward: 1.69
00300     //    Backward: 1.28
00301     //    Backward (long dt): 0.02
00302     //  Generic DoNewton:
00303     //    Forward: 1.66, 1.67, 1.67
00304     //    Backward: 1.36, 1.31, 1.32
00305     //    Backward (long dt): 0.02, 0.03, 0.03
00306     //  Generic DoNewton with JonW norm:
00307     //    Forward: 1.67, 1.67, 1.67
00308     //    Backward: 1.30, 1.29, 1.30
00309     //    Backward (long dt): 0.03, 0.03, 0.02
00310     //  CardiacNewtonSolver with JonW norm:
00311     //    Forward: 1.70, 1.69, 1.70, 1.71
00312     //    Backward: 1.29, 1.29, 1.30, 1.32
00313     //    Backward (long dt): 0.03, 0.03, 0.02, 0.03
00314     // So the method call itself has no significant impact, nor the use of generic solver.
00315     // The choice of norm also seems immaterial for this problem.
00316 }
00317 
00318 void BackwardEulerLuoRudyIModel1991::ComputeResidual(double var_environment__time, const double rCurrentGuess[1], double rResidual[1])
00319 {
00320     std::vector<double>& rY = rGetStateVariables();
00321     double membrane_V = rY[4];
00322     double slow_inward_current_d_gate_d = rY[5];
00323     double slow_inward_current_f_gate_f = rY[6];
00324 
00325     double slow_inward_current_E_si = 7.7-13.0287*log(rCurrentGuess[0]);
00326     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);
00327 
00328     double dCdt_using_current_guess = -1e-4*slow_inward_current_i_si+0.07*(1e-4-rCurrentGuess[0]);
00329 
00330     rResidual[0] = rCurrentGuess[0] - rY[3] - mDt*dCdt_using_current_guess;
00331 }
00332 
00333 void BackwardEulerLuoRudyIModel1991::ComputeJacobian(double var_environment__time, const double rCurrentGuess[1], double rJacobian[1][1])
00334 {
00335     std::vector<double>& rY = rGetStateVariables();
00336     double slow_inward_current_d_gate_d = rY[5];
00337     double slow_inward_current_f_gate_f = rY[6];
00338 
00339     rJacobian[0][0] = 1 - mDt*(-0.1172583e-3*slow_inward_current_d_gate_d*slow_inward_current_f_gate_f/rCurrentGuess[0]
00340                                - 0.07);
00341 }
00342 
00343 double BackwardEulerLuoRudyIModel1991::GetIIonic()
00344 {
00345     double fast_sodium_current_h_gate_h = mStateVariables[0];
00346     double fast_sodium_current_j_gate_j = mStateVariables[1];
00347     double fast_sodium_current_m_gate_m = mStateVariables[2];
00348     double intracellular_calcium_concentration_Cai = mStateVariables[3];
00349     double membrane_V = mStateVariables[4];
00350     double slow_inward_current_d_gate_d = mStateVariables[5];
00351     double slow_inward_current_f_gate_f = mStateVariables[6];
00352     double time_dependent_potassium_current_X_gate_X = mStateVariables[7];
00353 
00354     /*
00355      * Compute the LuoRudyIModel1991OdeSystem model
00356      */
00357     double background_current_i_b = background_current_g_b*(membrane_V-background_current_E_b);
00358 
00359     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);
00360     double slow_inward_current_E_si = 7.7-13.0287*log(intracellular_calcium_concentration_Cai);
00361 
00362     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);
00363 
00364     double time_dependent_potassium_current_g_K = 0.282*sqrt(ionic_concentrations_Ko/5.4);
00365     double time_dependent_potassium_current_Xi_gate_Xi;
00366 
00367     // Although the equation below looks strange (particularly the arguments of the
00368     // exponentials, it is in fact correct.
00369     if (membrane_V > -100.0)
00370     {
00371         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)));
00372     }
00373     else
00374     {
00375         #define COVERAGE_IGNORE
00376         time_dependent_potassium_current_Xi_gate_Xi = 1.0;
00377         #undef COVERAGE_IGNORE
00378     }
00379 
00380     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));
00381     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);
00382 
00383     double time_independent_potassium_current_g_K1 = 0.6047*sqrt(ionic_concentrations_Ko/5.4);
00384     double time_independent_potassium_current_E_K1 =((membrane_R*membrane_T)/membrane_F)*log(ionic_concentrations_Ko/ionic_concentrations_Ki);
00385     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)));
00386     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)));
00387     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);
00388     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);
00389 
00390     double plateau_potassium_current_Kp = 1.0/(1.0+exp((7.488-membrane_V)/5.98));
00391     double plateau_potassium_current_E_Kp = time_independent_potassium_current_E_K1;
00392     double plateau_potassium_current_i_Kp = plateau_potassium_current_g_Kp*plateau_potassium_current_Kp*(membrane_V-plateau_potassium_current_E_Kp);
00393 
00394     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;
00395 
00396     assert( !std::isnan(i_ionic));
00397     return i_ionic;
00398 }
00399 
00400 
00401 void BackwardEulerLuoRudyIModel1991::VerifyStateVariables()
00402 {
00403     const std::vector<double>& rY = rGetStateVariables();
00404 
00405     const double fast_sodium_current_h_gate_h = rY[0];            // gating
00406     const double fast_sodium_current_j_gate_j = rY[1];            // gating
00407     const double fast_sodium_current_m_gate_m = rY[2];            // gating
00408     const double intracellular_calcium_concentration_Cai = rY[3]; // concentration
00409     const double slow_inward_current_d_gate_d = rY[5];            // gating
00410     const double slow_inward_current_f_gate_f = rY[6];            // gating
00411     const double time_dependent_potassium_current_X_gate_X = rY[7]; // gating
00412 
00413     #define COVERAGE_IGNORE
00414     if (!(0.0<=fast_sodium_current_h_gate_h && fast_sodium_current_h_gate_h<=1.0))
00415     {
00416         EXCEPTION("h gate for fast sodium current has gone out of range. Check model parameters, for example spatial stepsize");
00417     }
00418 
00419     if (!(0.0<=fast_sodium_current_j_gate_j && fast_sodium_current_j_gate_j<=1.0))
00420     {
00421         EXCEPTION("j gate for fast sodium current has gone out of range. Check model parameters, for example spatial stepsize");
00422     }
00423 
00424     if (!(0.0<=fast_sodium_current_m_gate_m && fast_sodium_current_m_gate_m<=1.0))
00425     {
00426         EXCEPTION("m gate for fast sodium current has gone out of range. Check model parameters, for example spatial stepsize");
00427     }
00428 
00429     if (!(0.0<intracellular_calcium_concentration_Cai))
00430     {
00431         EXCEPTION("intracellular_calcium_concentration_Cai has become non-positive, ie gone out of range. Check model parameters, for example spatial stepsize");
00432     }
00433 
00434     if (!(0.0<=slow_inward_current_d_gate_d && slow_inward_current_d_gate_d<=1.0))
00435     {
00436         EXCEPTION("d gate for slow inward current has gone out of range. Check model parameters, for example spatial stepsize");
00437     }
00438 
00439     if (!(0.0<=slow_inward_current_f_gate_f && slow_inward_current_f_gate_f<=1.0))
00440     {
00441         EXCEPTION("f gate for slow inward current has gone out of range. Check model parameters, for example spatial stepsize");
00442     }
00443 
00444     if (!(0.0<=time_dependent_potassium_current_X_gate_X && time_dependent_potassium_current_X_gate_X<=1.0))
00445     {
00446         EXCEPTION("X gate for time dependent potassium current has gone out of range. Check model parameters, for example spatial stepsize");
00447     }
00448     #undef COVERAGE_IGNORE
00449 }
00450 
00451 
00452 
00453 // Serialization for Boost >= 1.36
00454 #include "SerializationExportWrapperForCpp.hpp"
00455 CHASTE_CLASS_EXPORT(BackwardEulerLuoRudyIModel1991)

Generated by  doxygen 1.6.2