DiFrancescoNoble1985OdeSystem.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 "DiFrancescoNoble1985OdeSystem.hpp"
00029 #include "OdeSystemInformation.hpp"
00030 #include <cmath>
00031 #include <cstdio>
00032 //#include <iostream>
00033 #include "Exception.hpp"
00034 
00035 //
00036 // Model-scope constant parameters
00037 //
00038 const double DiFrancescoNoble1985OdeSystem::membrane_C = 75.0;
00039 const double DiFrancescoNoble1985OdeSystem::radius = 0.05;
00040 const double DiFrancescoNoble1985OdeSystem::length = 2;
00041 const double DiFrancescoNoble1985OdeSystem::V_e_ratio = 0.1;
00042 const double DiFrancescoNoble1985OdeSystem::R = 8314.472;
00043 const double DiFrancescoNoble1985OdeSystem::T = 310;
00044 const double DiFrancescoNoble1985OdeSystem::F = 96485.3415;
00045 // conductances and the likes
00046 const double DiFrancescoNoble1985OdeSystem::g_fna = 3;
00047 const double DiFrancescoNoble1985OdeSystem::g_fk = 3;
00048 const double DiFrancescoNoble1985OdeSystem::g_na_b = 0.18;
00049 const double DiFrancescoNoble1985OdeSystem::g_ca_b = 0.02;
00050 const double DiFrancescoNoble1985OdeSystem::g_na = 750.0;
00051 const double DiFrancescoNoble1985OdeSystem::g_k1 = 920;
00052 const double DiFrancescoNoble1985OdeSystem::g_to = 0.28;
00053 const double DiFrancescoNoble1985OdeSystem::I_P = 125;
00054 const double DiFrancescoNoble1985OdeSystem::i_kmax = 180;
00055 const double DiFrancescoNoble1985OdeSystem::k_naca = 0.02;
00056 const double DiFrancescoNoble1985OdeSystem::P_si = 15.0;
00057 // concentrations
00058 const double DiFrancescoNoble1985OdeSystem::Nao = 140;
00059 const double DiFrancescoNoble1985OdeSystem::Cao = 2;
00060 const double DiFrancescoNoble1985OdeSystem::Kb = 4;
00061 const double DiFrancescoNoble1985OdeSystem::Kmf = 45;
00062 const double DiFrancescoNoble1985OdeSystem::Km1 = 210;
00063 const double DiFrancescoNoble1985OdeSystem::Kmto = 10;
00064 const double DiFrancescoNoble1985OdeSystem::KmK = 1;
00065 const double DiFrancescoNoble1985OdeSystem::KmNa = 40;
00066 const double DiFrancescoNoble1985OdeSystem::Kmf2 = 0.001;
00067 const double DiFrancescoNoble1985OdeSystem::KmCa = 0.001;
00068 const double DiFrancescoNoble1985OdeSystem::Km_Ca = 0.0005;
00069 // others
00070 const double DiFrancescoNoble1985OdeSystem::n_naca = 3;
00071 const double DiFrancescoNoble1985OdeSystem::gamma = 0.5;
00072 const double DiFrancescoNoble1985OdeSystem::d_naca = 0.001;
00073 const double DiFrancescoNoble1985OdeSystem::rCa = 2.0;
00074 const double DiFrancescoNoble1985OdeSystem::tau_up = 0.025;
00075 const double DiFrancescoNoble1985OdeSystem::tau_rep = 2;
00076 const double DiFrancescoNoble1985OdeSystem::tau_rel = 0.05;
00077 const double DiFrancescoNoble1985OdeSystem::Ca_up_max = 5;
00078 const double DiFrancescoNoble1985OdeSystem::pf = 0.0007;
00079 const double DiFrancescoNoble1985OdeSystem::Vecs = 0.05;
00080 
00081 
00082 /*Constructor*/
00083 DiFrancescoNoble1985OdeSystem::DiFrancescoNoble1985OdeSystem(
00084         boost::shared_ptr<AbstractIvpOdeSolver> pSolver,
00085         boost::shared_ptr<AbstractStimulusFunction> pIntracellularStimulus)
00086     : AbstractCardiacCell(pSolver, 16, 0, pIntracellularStimulus)
00087 {
00088     mpSystemInfo = OdeSystemInformation<DiFrancescoNoble1985OdeSystem>::Instance();
00089 
00090     //compute the fixed parameters
00091     Vcell= 3.141592654*pow(radius,2)*length;
00092     Vi=Vcell*(1.0-V_e_ratio);
00093     Vup=0.05*Vi;
00094     Vrel=0.02*Vi;
00095     Ve=Vi/(1.0-Vecs);
00096     RToNF=R*T/F;
00097 
00098     //initialise
00099     Init();
00100 }
00101 
00105 DiFrancescoNoble1985OdeSystem::~DiFrancescoNoble1985OdeSystem(void)
00106 {
00107 }
00108 
00118 void DiFrancescoNoble1985OdeSystem::EvaluateYDerivatives(double time,
00119                                                       const std::vector<double> &rY,
00120                                                       std::vector<double> &rDY)
00121 {
00122     double Vm = rY[0];
00123     double y_gate = rY[1];
00124     double s_gate = rY[2];
00125     double m_gate = rY[3];
00126     double h_gate = rY[4];
00127     double d_gate= rY[5];
00128     double f_gate = rY[6];
00129     double f2_gate = rY[7];
00130     double p = rY[8];
00131     double Nai = rY[9];
00132     double Cai = rY[10];
00133     double Kc = rY[11];
00134     double Ki = rY[12];
00135     double x_gate = rY[13];
00136     double Ca_up = rY[14];
00137     double Ca_rel = rY[15];
00138     VerifyStateVariables();
00139 
00141    // gating variables and ionic fluxes have been scaled by a factor of 1,000 to convert them into milliseconds
00143     double Emh=RToNF*log((Nao+0.12*Kc)/(Nai+0.12*Ki));
00144     double E_na=RToNF*log(Nao/Nai);
00145     double E_k=RToNF*log(Kc/Ki);
00146     double E_ca=0.5*RToNF*log(Cao/Cai);
00147     //i_f parameters
00148     double i_fk=(Kc/(Kc+Kmf))*(g_fk*(Vm-E_k));
00149     double i_fna=(Kc/(Kc+Kmf))*g_fna*(Vm-E_na);
00150     double alpha_y=0.001*0.05*exp(-0.067*(Vm+52.0-10.0));
00151     double beta_y;
00152 
00153    if(fabs(Vm+52.0-10.0)<=0.0001) //i.e. if Vm = (-52)
00154     {
00155         #define COVERAGE_IGNORE
00156         beta_y=2.5*0.001;
00157         #undef COVERAGE_IGNORE
00158     }
00159     else
00160     {
00161         beta_y=0.001*(Vm+52.0-10.0)/(1.0-exp(-0.2*(Vm+52.0-10.0)));
00162     }
00163     double y_gate_prime=alpha_y*(1.0-y_gate)-beta_y*y_gate;
00164 
00165     //i_k parameters
00166     double I_K=i_kmax*(Ki-Kc*exp(-Vm/RToNF))/140.0;
00167     double alpha_x=0.001*0.5*exp(0.0826*(Vm+50.0))/(1.0+exp(0.057*(Vm+50.0)));
00168     double beta_x=0.001*1.3*exp(-0.06*(Vm+20.0))/(1.0+exp(-0.04*(Vm+20.0)));
00169     double x_gate_prime=alpha_x*(1.0-x_gate)-beta_x*x_gate;
00170 
00171     //i_to parameters
00172     double alpha_s=0.001*0.033*exp(-Vm/17.0);
00173     double beta_s=0.001*33.0/(1.0+exp(-(Vm+10.0)/8.0));
00174     double s_gate_prime=alpha_s*(1.0-s_gate)-beta_s*s_gate;
00175 
00176     // i_na parameters
00177     double alpha_m;
00178     if(fabs(Vm+41.0)<=0.00001) //i.e. if Vm = -41
00179     {
00180       #define COVERAGE_IGNORE
00181       alpha_m=2000*0.001;
00182       #undef COVERAGE_IGNORE
00183     }
00184     else
00185     {
00186         alpha_m=0.001*200.0*(Vm+41.0)/(1.0-exp(-0.1*(Vm+41.0)));
00187     }
00188 
00189     double beta_m=0.001*8000*exp(-0.056*(Vm+66.0));
00190 
00191     double m_gate_prime=alpha_m*(1.0-m_gate)-beta_m*m_gate;
00192 
00193     double alpha_h=0.001*20*exp(-0.125*(Vm+75.0));
00194     double beta_h=0.001*2000.0/(320.0*exp(-0.1*(Vm+75.0))+1.0);
00195     double h_gate_prime=alpha_h*(1.0-h_gate)-beta_h*h_gate;
00196 
00197     //d gate
00198     double alpha_d;
00199     if(fabs(Vm+24.0-5.0)<=0.0001) //i.e. if Vm = -24
00200     {
00201         #define COVERAGE_IGNORE
00202         alpha_d=120*0.001;
00203         #undef COVERAGE_IGNORE
00204     }
00205     else
00206     {
00207         alpha_d=0.001*30*(Vm+24.0-5.0)/(1.0-exp(-(Vm+24.0-5.0)/4.0));
00208     }
00209 
00210     double beta_d;
00211     if(fabs(Vm+24.0-5.0)<=0.0001) //i.e. if Vm = -24
00212     {
00213         #define COVERAGE_IGNORE
00214         beta_d=120*0.001;
00215         #undef COVERAGE_IGNORE
00216     }
00217     else
00218     {
00219         beta_d=0.001*12.0*(Vm+24.0-5.0)/(exp((Vm+24.0-5.0)/10.0)-1.0);
00220     }
00221     double d_gate_prime=alpha_d*(1.0-d_gate)-beta_d*d_gate;
00222 
00223     //f gate.
00224     double alpha_f;
00225     if(fabs(Vm+34)<=0.0001)
00226     {
00227         #define COVERAGE_IGNORE
00228         alpha_f=25*0.001;
00229         #undef COVERAGE_IGNORE
00230     }
00231     else
00232     {
00233         alpha_f=0.001*6.25*(Vm+34.0)/(exp((Vm+34.0)/4.0)-1.0);
00234     }
00235     double beta_f=0.001*50.0/(1.0+exp(-(Vm+34.0)/4.0));
00236 
00237     double f_gate_prime=alpha_f*(1.0-f_gate)-beta_f*f_gate;
00238 
00239     //f2_gate
00240     double alpha_f2=0.001*5.0;
00241     double beta_f2=0.001*Cai*alpha_f2/Kmf2;
00242     double f2_gate_prime=alpha_f2-(f2_gate*(alpha_f2+beta_f2));
00243 
00244     //p_gate
00245     double alpha_p=0.001*(0.625*(Vm+34.0))/(exp((Vm+34.0)/4.0)-1.0);
00246     double beta_p=0.001*5.0/(1.0+exp(-(Vm+34.0)/4.0));
00247     double p_prime=alpha_p*(1.0-p)-beta_p*p;
00248 
00249 
00250       //i_na
00251     double i_na=g_na*pow(m_gate,3)*h_gate*(Vm-Emh);
00252 
00253     //i_f
00254     double i_f=y_gate*(i_fk+i_fna);
00255     //i_k
00256     double i_k=x_gate*I_K;
00257     //i_k1
00258     double i_k1=(g_k1*Kc/(Kc+Km1))*(Vm-E_k)/(1.0+exp((Vm-E_k+10)*2.0/RToNF));
00259     //i_to
00260     double i_to=((((s_gate*g_to*Cai*(0.2+Kc/(Kmto+Kc)))/(Km_Ca+Cai))*(Vm+10.0))/(1.0-exp(-0.2*(Vm+10.0))))*((Ki*exp((Vm*0.5)/RToNF))-(Kc*exp(-0.5*Vm/RToNF)));
00261     //i_na_b
00262     double i_na_b=g_na_b*(Vm-E_na);
00263     //i_p
00264     double i_p=((I_P*Kc)/(KmK+Kc))*Nai/(Nai+KmNa);
00265     //i_naca
00266     double i_naca=(k_naca*((Cao*(exp((gamma*Vm*(n_naca-2))/(RToNF)))*(pow(Nai,n_naca)))-(Cai*(exp(((gamma-1.0)*(n_naca-2.0)*Vm)/(RToNF)))*(pow(Nao,n_naca)))))/((1.0+(d_naca*((Cai*pow(Nao,n_naca))+(Cao*pow(Nai,n_naca)))))*(1.0+(Cai/0.0069)));
00267     //i_b_ca
00268     double i_ca_b=g_ca_b*(Vm-E_ca);
00269 
00270     //i_si
00271     double i_sica=((4.0*P_si*(Vm-50.0))/(RToNF*(1.0-(exp((-2.0*(Vm-50.0))/RToNF)))))*((Cai*exp(100.0/RToNF))-(Cao*exp(-2.0*(Vm-50.0)/RToNF)))*d_gate*f_gate*f2_gate;
00272     double i_sik=((0.01*P_si*(Vm-50.0))/(RToNF*(1.0-(exp((-(Vm-50.0))/RToNF)))))*((Ki*exp(50.0/RToNF))-(Kc*exp(-(Vm-50.0)/RToNF)))*d_gate*f_gate*f2_gate;
00273     double i_sina=((0.01*P_si*(Vm-50.0))/(RToNF*(1.0-(exp((-(Vm-50.0))/RToNF)))))*((Nai*exp(50.0/RToNF))-(Nao*exp(-(Vm-50.0)/RToNF)))*d_gate*f_gate*f2_gate;
00274     double i_si=i_sik+i_sica+i_sina;
00275 
00276     //Na ionic concentrations
00277     double Nai_prime= -0.001*(i_na+i_na_b+i_fna+i_sina+3*i_p+(n_naca*i_naca)/(n_naca-2.0))/(Vi*F);//mM/ms
00278 
00279     //Ca ionic concentrations
00280     double i_up=(2*F*Vi/(tau_up*Ca_up_max))*Cai*(Ca_up_max-Ca_up);//nA
00281     double i_tr=(2*F*Vrel/tau_rep)*p*(Ca_up-Ca_rel);//nA
00282     double i_rel=(2*F*Vrel/tau_rel)*Ca_rel*(pow(Cai,rCa))/(pow(Cai,rCa)+pow(KmCa,rCa));//nA
00283 
00284     double Ca_up_prime=0.001*(i_up-i_tr)/(2*Vup*F);//mM/ms
00285     double Ca_rel_prime=0.001*(i_tr-i_rel)/(2*Vrel*F);//mM/ms
00286     double Cai_prime=-0.001*(i_sica+i_ca_b+i_up-i_rel-2*(i_naca/(n_naca-2.0)))/(2*Vi*F);//mM/ms
00287 
00288     //K ionic concentrations
00289     double i_mk=(i_k1+i_k+i_fk+i_sik+i_to-2*i_p);//nA
00290     double Kc_prime=i_mk*0.001/(F*Ve)-pf*(Kc-Kb);//mM/ms
00291     double Ki_prime=(-i_mk*0.001/(Vi*F));//mM/ms
00292 
00293     //stimulus current
00294     double i_stim = GetStimulus(time);
00295 
00296     //calculate dV
00297     double Vm_prime = (-1.0/membrane_C)*(i_f+
00298                                          i_k+
00299                                          i_k1+
00300                                          i_to+
00301                                          i_na_b+
00302                                          i_p+
00303                                          i_naca+
00304                                          i_ca_b+
00305                                          i_na+
00306                                          i_si+
00307                                          i_stim);//mV/ms
00308     // do not update voltage if the mSetVoltageDerivativeToZero flag has been set
00309     if (mSetVoltageDerivativeToZero)
00310     {
00311         Vm_prime = 0;
00312     }
00313 
00314     rDY[0] = Vm_prime;
00315     rDY[1] = y_gate_prime;
00316     rDY[2] = s_gate_prime;
00317     rDY[3] = m_gate_prime;
00318     rDY[4] = h_gate_prime;
00319     rDY[5] = d_gate_prime;
00320     rDY[6] = f_gate_prime;
00321     rDY[7] = f2_gate_prime;
00322     rDY[8] = p_prime;
00323     rDY[9] = Nai_prime;
00324     rDY[10]= Cai_prime;
00325     rDY[11] = Kc_prime;
00326     rDY[12] = Ki_prime;
00327     rDY[13] = x_gate_prime;
00328     rDY[14] = Ca_up_prime;
00329     rDY[15] = Ca_rel_prime;
00330 }
00331 
00332 
00333 double DiFrancescoNoble1985OdeSystem::GetIIonic()
00334 {
00335 
00336     double Vm = mStateVariables[0];
00337     double y_gate = mStateVariables[1];
00338     double s_gate = mStateVariables[2];
00339     double m_gate = mStateVariables[3];
00340     double h_gate = mStateVariables[4];
00341     double d_gate= mStateVariables[5];
00342     double f_gate = mStateVariables[6];
00343     double f2_gate = mStateVariables[7];
00344     //double p = mStateVariables[8];
00345     double Nai = mStateVariables[9];
00346     double Cai = mStateVariables[10];
00347     double Kc = mStateVariables[11];
00348     double Ki = mStateVariables[12];
00349     double x_gate = mStateVariables[13];
00350     //double Ca_up = mStateVariables[14];
00351     //double Ca_rel = mStateVariables[15];
00352 
00353     /*
00354      * Compute the DiFrancescoNoble1985OdeSystem model
00355      */
00356     double Emh=RToNF*log((Nao+0.12*Kc)/(Nai+0.12*Ki));
00357     double E_na=RToNF*log(Nao/Nai);
00358     double E_k=RToNF*log(Kc/Ki);
00359     double E_ca=0.5*RToNF*log(Cao/Cai);
00360 
00361     //i_na
00362     double i_na=g_na*pow(m_gate,3)*h_gate*(Vm-Emh);
00363     //i_f
00364     double i_fk=(Kc/(Kc+Kmf))*(g_fk*(Vm-E_k));
00365     double i_fna=(Kc/(Kc+Kmf))*g_fna*(Vm-E_na);
00366     double i_f=y_gate*(i_fk+i_fna);
00367     //i_k
00368     double I_K=i_kmax*(Ki-Kc*exp(-Vm/RToNF))/140.0;
00369     double i_k=x_gate*I_K;
00370     //i_k1
00371     double i_k1=(g_k1*Kc/(Kc+Km1))*(Vm-E_k)/(1.0+exp((Vm-E_k+10)*2.0/RToNF));
00372     //i_to
00373     double i_to=((((s_gate*g_to*Cai*(0.2+Kc/(Kmto+Kc)))/(Km_Ca+Cai))*(Vm+10.0))/(1.0-exp(-0.2*(Vm+10.0))))*((Ki*exp((Vm*0.5)/RToNF))-(Kc*exp(-0.5*Vm/RToNF)));
00374     //i_na_b
00375     double i_na_b=g_na_b*(Vm-E_na);
00376     //i_p
00377     double i_p=((I_P*Kc)/(KmK+Kc))*Nai/(Nai+KmNa);
00378     //i_naca
00379     double i_naca=(k_naca*((Cao*(exp((gamma*Vm*(n_naca-2))/(RToNF)))*(pow(Nai,n_naca)))-(Cai*(exp(((gamma-1.0)*(n_naca-2.0)*Vm)/(RToNF)))*(pow(Nao,n_naca)))))/((1.0+(d_naca*((Cai*pow(Nao,n_naca))+(Cao*pow(Nai,n_naca)))))*(1.0+(Cai/0.0069)));
00380     //i_b_ca
00381     double i_ca_b=g_ca_b*(Vm-E_ca);
00382     //i_si
00383     double i_sica=((4.0*P_si*(Vm-50.0))/(RToNF*(1.0-(exp((-2.0*(Vm-50.0))/RToNF)))))*((Cai*exp(100.0/RToNF))-(Cao*exp(-2.0*(Vm-50.0)/RToNF)))*d_gate*f_gate*f2_gate;
00384     double i_sik=((0.01*P_si*(Vm-50.0))/(RToNF*(1.0-(exp((-(Vm-50.0))/RToNF)))))*((Ki*exp(50.0/RToNF))-(Kc*exp(-(Vm-50.0)/RToNF)))*d_gate*f_gate*f2_gate;
00385     double i_sina=((0.01*P_si*(Vm-50.0))/(RToNF*(1.0-(exp((-(Vm-50.0))/RToNF)))))*((Nai*exp(50.0/RToNF))-(Nao*exp(-(Vm-50.0)/RToNF)))*d_gate*f_gate*f2_gate;
00386     double i_si=i_sik+i_sica+i_sina;
00387 
00388     double i_ionic =   (i_f+
00389                          i_k+
00390                          i_k1+
00391                          i_to+
00392                          i_na_b+
00393                          i_p+
00394                          i_naca+
00395                          i_ca_b+
00396                          i_na+
00397                          i_si); /*this is in nA*/
00398 
00399     assert(!std::isnan(i_ionic));
00400 
00401     double i_ionic_in_microA_per_cm2=i_ionic*pow(10,-3)/0.075;
00402     return i_ionic_in_microA_per_cm2;
00403      /*   I_ion is in nA. I_ion*pow(10,-3) is in microA.
00404      *    Please note that in the mono/bidomain formulation, I_ion needs to be in microA/cm2.
00405      *    The cell capacitance is, from the paper, 0.075 microF.
00406      *    In order to get the cell density factor, we use the same method as in the Noble98 model and we obtain using the
00407      *    Cm=0.075microF from the cell model and Cm=1.0microF/cm^2 from the mono/bidomain equations.
00408      *    Hence, the cell density factor is (1.0 microF/cm2)/(0.075 microF).
00409      */
00410 }
00411 
00412 void DiFrancescoNoble1985OdeSystem::VerifyStateVariables()
00413 {
00414 //#ifndef NDEBUG
00415     const std::vector<double>& rY = rGetStateVariables();
00416 
00417     const double Vm = rY[0];
00418     const double y_gate = rY[1];
00419     const double r_gate = rY[2];
00420     const double m_gate = rY[3];
00421     const double h_gate = rY[4];
00422     const double d_gate= rY[5];
00423     const double f_gate = rY[6];
00424     const double f2_gate = rY[7];
00425     const double p = rY[8];
00426     const double Nai = rY[9];
00427     const double Cai = rY[10];
00428     const double Kc = rY[11];
00429     const double Ki = rY[12];
00430     const double x_gate = rY[13];
00431     const double Ca_up = rY[14];
00432     const double Ca_rel = rY[15];
00433     #define COVERAGE_IGNORE
00434     if (!(200>=Vm && Vm>=-200))
00435     {
00436         EXCEPTION(DumpState("Vm is really out of range!"));
00437     }
00438     if (!(0.0<=y_gate && y_gate<=1.0))
00439     {
00440         EXCEPTION(DumpState("y gate has gone out of range. Check model parameters, for example spatial stepsize"));
00441     }
00442     if (!(0.0<=r_gate && r_gate<=1.0))
00443     {
00444         EXCEPTION(DumpState("r gate has gone out of range. Check model parameters, for example spatial stepsize"));
00445     }
00446     if (!(0.0<=m_gate && m_gate<=1.0))
00447     {
00448         EXCEPTION(DumpState("m gate  has gone out of range. Check model parameters, for example spatial stepsize"));
00449     }
00450     if (!(0.0<=h_gate && h_gate<=1.0))
00451     {
00452         EXCEPTION(DumpState("h gate  has gone out of range. Check model parameters, for example spatial stepsize"));
00453     }
00454     if (!(0.0<=d_gate && d_gate<=1.0))
00455     {
00456         EXCEPTION(DumpState("d gate  has gone out of range. Check model parameters, for example spatial stepsize"));
00457     }
00458     if (!(0.0<=f_gate && f_gate<=1.0))
00459     {
00460         EXCEPTION(DumpState("f gate for the plateau current has gone out of range. Check model parameters, for example spatial stepsize"));
00461     }
00462     if (!(0.0<=f2_gate && f2_gate<=1.0))
00463     {
00464         EXCEPTION(DumpState("f2 gate has gone out of range. Check model parameters, for example spatial stepsize"));
00465     }
00466     if (!(0.0<=p && p<=1.0))
00467     {
00468         EXCEPTION(DumpState("p gate has gone out of range. Check model parameters, for example spatial stepsize"));
00469     }
00470     if (!(Nai>0.0))
00471     {
00472         EXCEPTION(DumpState("intarcellular Na concentration is negative!"));
00473     }
00474     if (!(Cai>0.0))
00475     {
00476         EXCEPTION(DumpState("intarcellular Ca concentration is negative!"));
00477     }
00478     if (!(Kc>0.0))
00479     {
00480         EXCEPTION(DumpState("extracellular K concentration is negative!"));
00481     }
00482     if (!(Ki>0.0))
00483     {
00484         EXCEPTION(DumpState("intarcellular K concentration is negative!"));
00485     }
00486     if (!(0.0<=x_gate && x_gate<=1.0))
00487     {
00488         EXCEPTION(DumpState("x gate  has gone out of range. Check model parameters, for example spatial stepsize"));
00489     }
00490     if (!(Ca_up>0.0))
00491     {
00492         EXCEPTION(DumpState("Ca_up concentration is negative!"));
00493     }
00494     if (!(Ca_rel>0.0))
00495     {
00496         EXCEPTION(DumpState("Ca_rel concentration is negative!"));
00497     }
00498     #undef COVERAGE_IGNORE
00499 //#endif
00500 }
00501 
00502 template<>
00503 void OdeSystemInformation<DiFrancescoNoble1985OdeSystem>::Initialise(void)
00504 {
00505     this->mVariableNames.push_back("V");
00506     this->mVariableUnits.push_back("mV");
00507     this->mInitialConditions.push_back(-87.01);
00508 
00509     this->mVariableNames.push_back("y_gate");
00510     this->mVariableUnits.push_back("");
00511     this->mInitialConditions.push_back(0.2);
00512 
00513     this->mVariableNames.push_back("s_gate");
00514     this->mVariableUnits.push_back("");
00515     this->mInitialConditions.push_back(1.0);
00516 
00517     this->mVariableNames.push_back("m_gate");
00518     this->mVariableUnits.push_back("");
00519     this->mInitialConditions.push_back(0.01);
00520 
00521     this->mVariableNames.push_back("h_gate");
00522     this->mVariableUnits.push_back("");
00523     this->mInitialConditions.push_back(0.8);
00524 
00525     this->mVariableNames.push_back("d_gate");
00526     this->mVariableUnits.push_back("");
00527     this->mInitialConditions.push_back(0.005);
00528 
00529     this->mVariableNames.push_back("f_gate");
00530     this->mVariableUnits.push_back("");
00531     this->mInitialConditions.push_back(1.0);
00532 
00533     this->mVariableNames.push_back("f2_gate");
00534     this->mVariableUnits.push_back("");
00535     this->mInitialConditions.push_back(1.0);
00536 
00537     this->mVariableNames.push_back("p");
00538     this->mVariableUnits.push_back("");
00539     this->mInitialConditions.push_back(1.0);
00540 
00541     this->mVariableNames.push_back("Nai");
00542     this->mVariableUnits.push_back("");
00543     this->mInitialConditions.push_back(8.0);
00544 
00545     this->mVariableNames.push_back("Cai");
00546     this->mVariableUnits.push_back("");
00547     this->mInitialConditions.push_back(0.00005);
00548 
00549     this->mVariableNames.push_back("Kc");
00550     this->mVariableUnits.push_back("");
00551     this->mInitialConditions.push_back(4.0);
00552 
00553     this->mVariableNames.push_back("Ki");
00554     this->mVariableUnits.push_back("");
00555     this->mInitialConditions.push_back(140.0);
00556 
00557     this->mVariableNames.push_back("x_gate");
00558     this->mVariableUnits.push_back("");
00559     this->mInitialConditions.push_back(0.01);
00560 
00561     this->mVariableNames.push_back("Ca_up");
00562     this->mVariableUnits.push_back("");
00563     this->mInitialConditions.push_back(2.0);
00564 
00565     this->mVariableNames.push_back("Ca_rel");
00566     this->mVariableUnits.push_back("");
00567     this->mInitialConditions.push_back(1.0);
00568 
00569     this->mInitialised = true;
00570 }

Generated by  doxygen 1.6.2