TenTusscher2006OdeSystem.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 "TenTusscher2006OdeSystem.hpp"
00029 #include "OdeSystemInformation.hpp"
00030 #include <cmath>
00031 #include <cstdio>
00032 #include "Exception.hpp"
00033 
00034 
00035 //
00036 // Model-scope constant parameters
00037 //
00038 const double TenTusscher2006OdeSystem::L_type_Ca_current_g_CaL = 0.0000398;   // nanoS_per_picoF
00039 const double TenTusscher2006OdeSystem::calcium_background_current_g_bca = 0.000592;   // nanoS_per_picoF
00040 const double TenTusscher2006OdeSystem::calcium_dynamics_Buf_c = 0.2;   // millimolar
00041 const double TenTusscher2006OdeSystem::calcium_dynamics_Buf_sr = 10.0;   // millimolar
00042 const double TenTusscher2006OdeSystem::calcium_dynamics_Buf_ss = 0.4;   // millimolar
00043 const double TenTusscher2006OdeSystem::calcium_dynamics_Ca_o = 2.0;   // millimolar
00044 const double TenTusscher2006OdeSystem::calcium_dynamics_EC = 1.5;   // millimolar
00045 const double TenTusscher2006OdeSystem::calcium_dynamics_K_buf_c = 0.001;   // millimolar
00046 const double TenTusscher2006OdeSystem::calcium_dynamics_K_buf_sr = 0.3;   // millimolar
00047 const double TenTusscher2006OdeSystem::calcium_dynamics_K_buf_ss = 0.00025;   // millimolar
00048 const double TenTusscher2006OdeSystem::calcium_dynamics_K_up = 0.00025;   // millimolar
00049 const double TenTusscher2006OdeSystem::calcium_dynamics_V_leak = 0.00036;   // millimolar_per_millisecond
00050 const double TenTusscher2006OdeSystem::calcium_dynamics_V_rel = 0.102;   // millimolar_per_millisecond
00051 const double TenTusscher2006OdeSystem::calcium_dynamics_V_sr = 0.001094;   // micrometre3
00052 const double TenTusscher2006OdeSystem::calcium_dynamics_V_ss = 0.00005468;   // micrometre3
00053 const double TenTusscher2006OdeSystem::calcium_dynamics_V_xfer = 0.0038;   // millimolar_per_millisecond
00054 const double TenTusscher2006OdeSystem::calcium_dynamics_Vmax_up = 0.006375;   // millimolar_per_millisecond
00055 const double TenTusscher2006OdeSystem::calcium_dynamics_k1_prime = 0.15;   // per_millimolar2_per_millisecond
00056 const double TenTusscher2006OdeSystem::calcium_dynamics_k2_prime = 0.045;   // per_millimolar_per_millisecond
00057 const double TenTusscher2006OdeSystem::calcium_dynamics_k3 = 0.06;   // per_millisecond
00058 const double TenTusscher2006OdeSystem::calcium_dynamics_k4 = 0.005;   // per_millisecond
00059 const double TenTusscher2006OdeSystem::calcium_dynamics_max_sr = 2.5;   // dimensionless
00060 const double TenTusscher2006OdeSystem::calcium_dynamics_min_sr = 1.0;   // dimensionless
00061 const double TenTusscher2006OdeSystem::calcium_pump_current_K_pCa = 0.0005;   // millimolar
00062 const double TenTusscher2006OdeSystem::calcium_pump_current_g_pCa = 0.1238;   // nanoS_per_picoF
00063 const double TenTusscher2006OdeSystem::fast_sodium_current_g_Na = 14.838;   // nanoS_per_picoF
00064 const double TenTusscher2006OdeSystem::inward_rectifier_potassium_current_g_K1 = 5.405;   // nanoS_per_picoF
00065 const double TenTusscher2006OdeSystem::membrane_Cm = 0.185;   // microF_per_cm2
00066 const double TenTusscher2006OdeSystem::membrane_F = 96485.3415;   // coulomb_per_millimole
00067 const double TenTusscher2006OdeSystem::membrane_R = 8314.472;   // joule_per_mole_kelvin
00068 const double TenTusscher2006OdeSystem::membrane_T = 310.0;   // kelvin
00069 const double TenTusscher2006OdeSystem::membrane_V_c = 0.016404;   // micrometre3
00070 const double TenTusscher2006OdeSystem::potassium_dynamics_K_o = 5.4;   // millimolar
00071 const double TenTusscher2006OdeSystem::potassium_pump_current_g_pK = 0.0146;   // nanoS_per_picoF
00072 const double TenTusscher2006OdeSystem::rapid_time_dependent_potassium_current_g_Kr = 0.153;   // nanoS_per_picoF
00073 const double TenTusscher2006OdeSystem::reversal_potentials_P_kna = 0.03;   // nanoA_per_millimolar
00074 const double TenTusscher2006OdeSystem::slow_time_dependent_potassium_current_g_Ks = 0.392;   // nanoS_per_picoF
00075 const double TenTusscher2006OdeSystem::sodium_background_current_g_bna = 0.00029;   // nanoS_per_picoF
00076 const double TenTusscher2006OdeSystem::sodium_calcium_exchanger_current_K_NaCa = 1000.0;   // picoA_per_picoF
00077 const double TenTusscher2006OdeSystem::sodium_calcium_exchanger_current_K_sat = 0.1;   // dimensionless
00078 const double TenTusscher2006OdeSystem::sodium_calcium_exchanger_current_Km_Ca = 1.38;   // millimolar
00079 const double TenTusscher2006OdeSystem::sodium_calcium_exchanger_current_Km_Nai = 87.5;   // millimolar
00080 const double TenTusscher2006OdeSystem::sodium_calcium_exchanger_current_alpha = 2.5;   // dimensionless
00081 const double TenTusscher2006OdeSystem::sodium_calcium_exchanger_current_gamma = 0.35;   // dimensionless
00082 const double TenTusscher2006OdeSystem::sodium_dynamics_Na_o = 140.0;   // millimolar
00083 const double TenTusscher2006OdeSystem::sodium_potassium_pump_current_K_mNa = 40.0;   // millimolar
00084 const double TenTusscher2006OdeSystem::sodium_potassium_pump_current_K_mk = 1.0;   // millimolar
00085 const double TenTusscher2006OdeSystem::sodium_potassium_pump_current_P_NaK = 2.724;   // picoA_per_picoF
00086 const double TenTusscher2006OdeSystem::transient_outward_current_g_to = 0.294;   // nanoS_per_picoF
00087 
00088 
00089 
00090 /*Constructor*/
00091 TenTusscher2006OdeSystem::TenTusscher2006OdeSystem(
00092         boost::shared_ptr<AbstractIvpOdeSolver> pSolver,
00093         boost::shared_ptr<AbstractStimulusFunction> pIntracellularStimulus)
00094     : AbstractCardiacCell(pSolver, 19, 11, pIntracellularStimulus)
00095 {
00096 
00097     mpSystemInfo = OdeSystemInformation<TenTusscher2006OdeSystem>::Instance();
00098     //Initialise the scale factors
00099     mScaleFactorGks=1.0;
00100     mScaleFactorIto=1.0;
00101     mScaleFactorGkr=1.0;
00102 
00103     Init();
00104 }
00105 
00106 //Destructor
00107 TenTusscher2006OdeSystem::~TenTusscher2006OdeSystem(void)
00108 {
00109 }
00110 
00111 void TenTusscher2006OdeSystem::SetScaleFactorGks(double sfgks)
00112 {
00113     mScaleFactorGks=sfgks;
00114 }
00115 void TenTusscher2006OdeSystem::SetScaleFactorIto(double sfito)
00116 {
00117     mScaleFactorIto=sfito;
00118 }
00119 void TenTusscher2006OdeSystem::SetScaleFactorGkr(double sfgkr)
00120 {
00121     mScaleFactorGkr=sfgkr;
00122 }
00123 
00124 void TenTusscher2006OdeSystem::EvaluateYDerivatives(double time,
00125                                                     const std::vector<double> &rY,
00126                                                     std::vector<double> &rDY)
00127 {
00128    //---------------------------------------------------------------------------
00129    // State variables, initial value and names in comments beside
00130    //---------------------------------------------------------------------------
00131       //Vector for state variables//
00132       double Y[19];
00133       double dY[19];
00134 
00135     Y[0] = rY[0];// 3.373e-5;   // L_t Ype_Ca_current_d_gate_d (dimensionless)
00136     Y[1] = rY[1];// 0.9755;   // L_t Ype_Ca_current_f2_gate_f2 (dimensionless)
00137     Y[2] = rY[2];// 0.9953;   // L_t Ype_Ca_current_fCass_gate_fCass (dimensionless)
00138     Y[3] = rY[3];// 0.7888;   // L_t Ype_Ca_current_f_gate_f (dimensionless)
00139     Y[4] = rY[4];// 3.64;   // calcium_d Ynamics_Ca_SR (millimolar)
00140     Y[5] = rY[5];// 0.000126;   // calcium_d Ynamics_Ca_i (millimolar)
00141     Y[6] = rY[6];// 0.00036;   // calcium_d Ynamics_Ca_ss (millimolar)
00142     Y[7] = rY[7];// 0.9073;   // calcium_d Ynamics_R_prime (dimensionless)
00143     Y[8] = rY[8];// 0.7444;   // fast_sodium_current_h_gate_h (dimensionless)
00144     Y[9] = rY[9];// 0.7045;   // fast_sodium_current_j_gate_j (dimensionless)
00145     Y[10] = rY[10];// 0.00172;   // fast_sodium_current_m_gate_m (dimensionless)
00146     Y[11] = rY[11];// -85.23;   // membrane_V (millivolt)
00147     Y[12] = rY[12];// 136.89;   // potassium_d Ynamics_K_i (millimolar)
00148     Y[13] = rY[13];// 0.00621;   // rapid_time_dependent_potassium_current_Xr1_gate_Xr1 (dimensionless)
00149     Y[14] = rY[14];// 0.4712;   // rapid_time_dependent_potassium_current_Xr2_gate_Xr2 (dimensionless)
00150     Y[15] = rY[15];// 0.0095;   // slow_time_dependent_potassium_current_Xs_gate_Xs (dimensionless)
00151     Y[16] = rY[16];// 8.604;   // sodium_d Ynamics_Na_i (millimolar)
00152     Y[17] = rY[17];// 2.42e-8;   // transient_outward_current_r_gate_r (dimensionless)
00153     Y[18] = rY[18];// 0.999998;   // transient_outward_current_s_gate_s (dimensionless)
00154 
00155     VerifyStateVariables();
00156 
00157    double L_type_Ca_current_i_CaL = L_type_Ca_current_g_CaL*Y[0]*Y[3]*Y[1]*Y[2]*4.0*(Y[11]-15.0)*pow(membrane_F, 2.0)/(membrane_R*membrane_T)*(0.25*Y[6]*exp(2.0*(Y[11]-15.0)*membrane_F/(membrane_R*membrane_T))-calcium_dynamics_Ca_o)/(exp(2.0*(Y[11]-15.0)*membrane_F/(membrane_R*membrane_T))-1.0);
00158    double L_type_Ca_current_d_gate_d_inf = 1.0/(1.0+exp((-8.0-Y[11])/7.5));
00159    double L_type_Ca_current_d_gate_alpha_d = 1.4/(1.0+exp((-35.0-Y[11])/13.0))+0.25;
00160    double L_type_Ca_current_d_gate_beta_d = 1.4/(1.0+exp((Y[11]+5.0)/5.0));
00161    double L_type_Ca_current_d_gate_gamma_d = 1.0/(1.0+exp((50.0-Y[11])/20.0));
00162    double L_type_Ca_current_d_gate_tau_d = L_type_Ca_current_d_gate_alpha_d*L_type_Ca_current_d_gate_beta_d+L_type_Ca_current_d_gate_gamma_d;
00163    dY[0] = (L_type_Ca_current_d_gate_d_inf-Y[0])/L_type_Ca_current_d_gate_tau_d;
00164    double L_type_Ca_current_f2_gate_f2_inf = 0.67/(1.0+exp((Y[11]+35.0)/7.0))+0.33;
00165    double L_type_Ca_current_f2_gate_tau_f2 = 562.0*exp(-pow(Y[11]+27.0, 2.0)/240.0)+31.0/(1.0+exp((25.0-Y[11])/10.0))+80.0/(1.0+exp((Y[11]+30.0)/10.0));
00166    dY[1] = (L_type_Ca_current_f2_gate_f2_inf-Y[1])/L_type_Ca_current_f2_gate_tau_f2;
00167    double L_type_Ca_current_fCass_gate_fCass_inf = 0.6/(1.0+pow(Y[6]/0.05, 2.0))+0.4;
00168    double L_type_Ca_current_fCass_gate_tau_fCass = 80.0/(1.0+pow(Y[6]/0.05, 2.0))+2.0;
00169    dY[2] = (L_type_Ca_current_fCass_gate_fCass_inf-Y[2])/L_type_Ca_current_fCass_gate_tau_fCass;
00170    double L_type_Ca_current_f_gate_f_inf = 1.0/(1.0+exp((Y[11]+20.0)/7.0));
00171    double L_type_Ca_current_f_gate_tau_f = 1102.5*exp(-pow(Y[11]+27.0, 2.0)/225.0)+200.0/(1.0+exp((13.0-Y[11])/10.0))+180.0/(1.0+exp((Y[11]+30.0)/10.0))+20.0;
00172    dY[3] = (L_type_Ca_current_f_gate_f_inf-Y[3])/L_type_Ca_current_f_gate_tau_f;
00173    double reversal_potentials_E_Ca = 0.5*membrane_R*membrane_T/membrane_F*log(calcium_dynamics_Ca_o/Y[5]);
00174    double calcium_background_current_i_b_Ca = calcium_background_current_g_bca*(Y[11]-reversal_potentials_E_Ca);
00175    double calcium_dynamics_kcasr = calcium_dynamics_max_sr-(calcium_dynamics_max_sr-calcium_dynamics_min_sr)/(1.0+pow(calcium_dynamics_EC/Y[4], 2.0));
00176    double calcium_dynamics_k1 = calcium_dynamics_k1_prime/calcium_dynamics_kcasr;
00177    double calcium_dynamics_O = calcium_dynamics_k1*pow(Y[6], 2.0)*Y[7]/(calcium_dynamics_k3+calcium_dynamics_k1*pow(Y[6], 2.0));
00178    double calcium_dynamics_i_rel = calcium_dynamics_V_rel*calcium_dynamics_O*(Y[4]-Y[6]);
00179    double calcium_dynamics_i_up = calcium_dynamics_Vmax_up/(1.0+pow(calcium_dynamics_K_up, 2.0)/pow(Y[5], 2.0));
00180    double calcium_dynamics_i_leak = calcium_dynamics_V_leak*(Y[4]-Y[5]);
00181    double calcium_dynamics_i_xfer = calcium_dynamics_V_xfer*(Y[6]-Y[5]);
00182    double calcium_dynamics_k2 = calcium_dynamics_k2_prime*calcium_dynamics_kcasr;
00183    dY[7] = -calcium_dynamics_k2*Y[6]*Y[7]+calcium_dynamics_k4*(1.0-Y[7]);
00184    double calcium_dynamics_Ca_i_bufc = 1.0/(1.0+calcium_dynamics_Buf_c*calcium_dynamics_K_buf_c/pow(Y[5]+calcium_dynamics_K_buf_c, 2.0));
00185    double calcium_dynamics_Ca_sr_bufsr = 1.0/(1.0+calcium_dynamics_Buf_sr*calcium_dynamics_K_buf_sr/pow(Y[4]+calcium_dynamics_K_buf_sr, 2.0));
00186    double calcium_dynamics_Ca_ss_bufss = 1.0/(1.0+calcium_dynamics_Buf_ss*calcium_dynamics_K_buf_ss/pow(Y[6]+calcium_dynamics_K_buf_ss, 2.0));
00187    double calcium_pump_current_i_p_Ca = calcium_pump_current_g_pCa*Y[5]/(Y[5]+calcium_pump_current_K_pCa);
00188    double sodium_calcium_exchanger_current_i_NaCa = sodium_calcium_exchanger_current_K_NaCa*(exp(sodium_calcium_exchanger_current_gamma*Y[11]*membrane_F/(membrane_R*membrane_T))*pow(Y[16], 3.0)*calcium_dynamics_Ca_o-exp((sodium_calcium_exchanger_current_gamma-1.0)*Y[11]*membrane_F/(membrane_R*membrane_T))*pow(sodium_dynamics_Na_o, 3.0)*Y[5]*sodium_calcium_exchanger_current_alpha)/((pow(sodium_calcium_exchanger_current_Km_Nai, 3.0)+pow(sodium_dynamics_Na_o, 3.0))*(sodium_calcium_exchanger_current_Km_Ca+calcium_dynamics_Ca_o)*(1.0+sodium_calcium_exchanger_current_K_sat*exp((sodium_calcium_exchanger_current_gamma-1.0)*Y[11]*membrane_F/(membrane_R*membrane_T))));
00189    dY[5] = calcium_dynamics_Ca_i_bufc*((calcium_dynamics_i_leak-calcium_dynamics_i_up)*calcium_dynamics_V_sr/membrane_V_c+calcium_dynamics_i_xfer-(calcium_background_current_i_b_Ca+calcium_pump_current_i_p_Ca-2.0*sodium_calcium_exchanger_current_i_NaCa)*membrane_Cm/(2.0*membrane_V_c*membrane_F));
00190    dY[4] = calcium_dynamics_Ca_sr_bufsr*(calcium_dynamics_i_up-(calcium_dynamics_i_rel+calcium_dynamics_i_leak));
00191    dY[6] = calcium_dynamics_Ca_ss_bufss*(-L_type_Ca_current_i_CaL*membrane_Cm/(2.0*calcium_dynamics_V_ss*membrane_F)+calcium_dynamics_i_rel*calcium_dynamics_V_sr/calcium_dynamics_V_ss-calcium_dynamics_i_xfer*membrane_V_c/calcium_dynamics_V_ss);
00192    double reversal_potentials_E_Na = membrane_R*membrane_T/membrane_F*log(sodium_dynamics_Na_o/Y[16]);
00193    double fast_sodium_current_i_Na = fast_sodium_current_g_Na*pow(Y[10], 3.0)*Y[8]*Y[9]*(Y[11]-reversal_potentials_E_Na);
00194    double fast_sodium_current_h_gate_h_inf = 1.0/pow(1.0+exp((Y[11]+71.55)/7.43), 2.0);
00195 
00196    double fast_sodium_current_h_gate_alpha_h,fast_sodium_current_h_gate_beta_h;
00197    if (Y[11] < -40.0)
00198       fast_sodium_current_h_gate_alpha_h = 0.057*exp(-(Y[11]+80.0)/6.8);
00199    else
00200       fast_sodium_current_h_gate_alpha_h = 0.0;
00201 
00202    if (Y[11] < -40.0)
00203       fast_sodium_current_h_gate_beta_h = 2.7*exp(0.079*Y[11])+310000.0*exp(0.3485*Y[11]);
00204    else
00205       fast_sodium_current_h_gate_beta_h = 0.77/(0.13*(1.0+exp((Y[11]+10.66)/-11.1)));
00206 
00207    double fast_sodium_current_h_gate_tau_h = 1.0/(fast_sodium_current_h_gate_alpha_h+fast_sodium_current_h_gate_beta_h);
00208    dY[8] = (fast_sodium_current_h_gate_h_inf-Y[8])/fast_sodium_current_h_gate_tau_h;
00209    double fast_sodium_current_j_gate_j_inf = 1.0/pow(1.0+exp((Y[11]+71.55)/7.43), 2.0);
00210 
00211    double fast_sodium_current_j_gate_alpha_j ,fast_sodium_current_j_gate_beta_j;
00212    if (Y[11] < -40.0)
00213       fast_sodium_current_j_gate_alpha_j = (-25428.0*exp(0.2444*Y[11])-6.948e-6*exp(-0.04391*Y[11]))*(Y[11]+37.78)/(1.0+exp(0.311*(Y[11]+79.23)));
00214    else
00215       fast_sodium_current_j_gate_alpha_j = 0.0;
00216 
00217    if (Y[11] < -40.0)
00218       fast_sodium_current_j_gate_beta_j = 0.02424*exp(-0.01052*Y[11])/(1.0+exp(-0.1378*(Y[11]+40.14)));
00219    else
00220       fast_sodium_current_j_gate_beta_j = 0.6*exp(0.057*Y[11])/(1.0+exp(-0.1*(Y[11]+32.0)));
00221 
00222    double fast_sodium_current_j_gate_tau_j = 1.0/(fast_sodium_current_j_gate_alpha_j+fast_sodium_current_j_gate_beta_j);
00223    dY[9] = (fast_sodium_current_j_gate_j_inf-Y[9])/fast_sodium_current_j_gate_tau_j;
00224    double fast_sodium_current_m_gate_m_inf = 1.0/pow(1.0+exp((-56.86-Y[11])/9.03), 2.0);
00225    double fast_sodium_current_m_gate_alpha_m = 1.0/(1.0+exp((-60.0-Y[11])/5.0));
00226    double fast_sodium_current_m_gate_beta_m = 0.1/(1.0+exp((Y[11]+35.0)/5.0))+0.1/(1.0+exp((Y[11]-50.0)/200.0));
00227    double fast_sodium_current_m_gate_tau_m = fast_sodium_current_m_gate_alpha_m*fast_sodium_current_m_gate_beta_m;
00228    dY[10] = (fast_sodium_current_m_gate_m_inf-Y[10])/fast_sodium_current_m_gate_tau_m;
00229    double reversal_potentials_E_K = membrane_R*membrane_T/membrane_F*log(potassium_dynamics_K_o/Y[12]);
00230    double inward_rectifier_potassium_current_alpha_K1 = 0.1/(1.0+exp(0.06*(Y[11]-reversal_potentials_E_K-200.0)));
00231    double inward_rectifier_potassium_current_beta_K1 = (3.0*exp(0.0002*(Y[11]-reversal_potentials_E_K+100.0))+exp(0.1*(Y[11]-reversal_potentials_E_K-10.0)))/(1.0+exp(-0.5*(Y[11]-reversal_potentials_E_K)));
00232    double inward_rectifier_potassium_current_xK1_inf = inward_rectifier_potassium_current_alpha_K1/(inward_rectifier_potassium_current_alpha_K1+inward_rectifier_potassium_current_beta_K1);
00233    double inward_rectifier_potassium_current_i_K1 = inward_rectifier_potassium_current_g_K1*inward_rectifier_potassium_current_xK1_inf*(Y[11]-reversal_potentials_E_K);
00234 
00235 
00236    double transient_outward_current_i_to = mScaleFactorIto*transient_outward_current_g_to*Y[17]*Y[18]*(Y[11]-reversal_potentials_E_K);
00237    double rapid_time_dependent_potassium_current_i_Kr = mScaleFactorGkr*rapid_time_dependent_potassium_current_g_Kr*sqrt(potassium_dynamics_K_o/5.4)*Y[13]*Y[14]*(Y[11]-reversal_potentials_E_K);
00238    double reversal_potentials_E_Ks = membrane_R*membrane_T/membrane_F*log((potassium_dynamics_K_o+reversal_potentials_P_kna*sodium_dynamics_Na_o)/(Y[12]+reversal_potentials_P_kna*Y[16]));
00239    double slow_time_dependent_potassium_current_i_Ks = mScaleFactorGks*slow_time_dependent_potassium_current_g_Ks*pow(Y[15], 2.0)*(Y[11]-reversal_potentials_E_Ks);
00240    double sodium_potassium_pump_current_i_NaK = sodium_potassium_pump_current_P_NaK*potassium_dynamics_K_o/(potassium_dynamics_K_o+sodium_potassium_pump_current_K_mk)*Y[16]/(Y[16]+sodium_potassium_pump_current_K_mNa)/(1.0+0.1245*exp(-0.1*Y[11]*membrane_F/(membrane_R*membrane_T))+0.0353*exp(-Y[11]*membrane_F/(membrane_R*membrane_T)));
00241    double sodium_background_current_i_b_Na = sodium_background_current_g_bna*(Y[11]-reversal_potentials_E_Na);
00242    double potassium_pump_current_i_p_K = potassium_pump_current_g_pK*(Y[11]-reversal_potentials_E_K)/(1.0+exp((25.0-Y[11])/5.98));
00243 
00244    //stimulus current
00245    double i_stim = GetStimulus(time);
00246 
00247    dY[11] = -1.0/1.0*(inward_rectifier_potassium_current_i_K1+transient_outward_current_i_to+rapid_time_dependent_potassium_current_i_Kr+slow_time_dependent_potassium_current_i_Ks+L_type_Ca_current_i_CaL+sodium_potassium_pump_current_i_NaK+fast_sodium_current_i_Na+sodium_background_current_i_b_Na+sodium_calcium_exchanger_current_i_NaCa+calcium_background_current_i_b_Ca+potassium_pump_current_i_p_K+calcium_pump_current_i_p_Ca+i_stim);
00248    dY[12] = -(inward_rectifier_potassium_current_i_K1+transient_outward_current_i_to+rapid_time_dependent_potassium_current_i_Kr+slow_time_dependent_potassium_current_i_Ks+potassium_pump_current_i_p_K+i_stim-2.0*sodium_potassium_pump_current_i_NaK)/(membrane_V_c*membrane_F)*membrane_Cm;
00249    double rapid_time_dependent_potassium_current_Xr1_gate_xr1_inf = 1.0/(1.0+exp((-26.0-Y[11])/7.0));
00250    double rapid_time_dependent_potassium_current_Xr1_gate_alpha_xr1 = 450.0/(1.0+exp((-45.0-Y[11])/10.0));
00251    double rapid_time_dependent_potassium_current_Xr1_gate_beta_xr1 = 6.0/(1.0+exp((Y[11]+30.0)/11.5));
00252    double rapid_time_dependent_potassium_current_Xr1_gate_tau_xr1 = rapid_time_dependent_potassium_current_Xr1_gate_alpha_xr1*rapid_time_dependent_potassium_current_Xr1_gate_beta_xr1;
00253    dY[13] = (rapid_time_dependent_potassium_current_Xr1_gate_xr1_inf-Y[13])/rapid_time_dependent_potassium_current_Xr1_gate_tau_xr1;
00254    double rapid_time_dependent_potassium_current_Xr2_gate_xr2_inf = 1.0/(1.0+exp((Y[11]+88.0)/24.0));
00255    double rapid_time_dependent_potassium_current_Xr2_gate_alpha_xr2 = 3.0/(1.0+exp((-60.0-Y[11])/20.0));
00256    double rapid_time_dependent_potassium_current_Xr2_gate_beta_xr2 = 1.12/(1.0+exp((Y[11]-60.0)/20.0));
00257    double rapid_time_dependent_potassium_current_Xr2_gate_tau_xr2 = rapid_time_dependent_potassium_current_Xr2_gate_alpha_xr2*rapid_time_dependent_potassium_current_Xr2_gate_beta_xr2;
00258    dY[14] = (rapid_time_dependent_potassium_current_Xr2_gate_xr2_inf-Y[14])/rapid_time_dependent_potassium_current_Xr2_gate_tau_xr2;
00259    double slow_time_dependent_potassium_current_Xs_gate_xs_inf = 1.0/(1.0+exp((-5.0-Y[11])/14.0));
00260    double slow_time_dependent_potassium_current_Xs_gate_alpha_xs = 1400.0/sqrt(1.0+exp((5.0-Y[11])/6.0));
00261    double slow_time_dependent_potassium_current_Xs_gate_beta_xs = 1.0/(1.0+exp((Y[11]-35.0)/15.0));
00262    double slow_time_dependent_potassium_current_Xs_gate_tau_xs = slow_time_dependent_potassium_current_Xs_gate_alpha_xs*slow_time_dependent_potassium_current_Xs_gate_beta_xs+80.0;
00263    dY[15] = (slow_time_dependent_potassium_current_Xs_gate_xs_inf-Y[15])/slow_time_dependent_potassium_current_Xs_gate_tau_xs;
00264    dY[16] = -(fast_sodium_current_i_Na+sodium_background_current_i_b_Na+3.0*sodium_potassium_pump_current_i_NaK+3.0*sodium_calcium_exchanger_current_i_NaCa)/(membrane_V_c*membrane_F)*membrane_Cm;
00265    double transient_outward_current_r_gate_r_inf = 1.0/(1.0+exp((20.0-Y[11])/6.0));
00266    double transient_outward_current_r_gate_tau_r = 9.5*exp(-pow(Y[11]+40.0, 2.0)/1800.0)+0.8;
00267    dY[17] = (transient_outward_current_r_gate_r_inf-Y[17])/transient_outward_current_r_gate_tau_r;
00268    double transient_outward_current_s_gate_s_inf = 1.0/(1.0+exp((Y[11]+20.0)/5.0));
00269    double transient_outward_current_s_gate_tau_s = 85.0*exp(-pow(Y[11]+45.0, 2.0)/320.0)+5.0/(1.0+exp((Y[11]-20.0)/5.0))+3.0;
00270    dY[18] = (transient_outward_current_s_gate_s_inf-Y[18])/transient_outward_current_s_gate_tau_s;
00271 
00272 
00273     // do not update voltage if the mSetVoltageDerivativeToZero flag has been set
00274     if (mSetVoltageDerivativeToZero)
00275     {
00276         dY[11] = 0;
00277     }
00278 
00279     rDY[0] = dY[0];
00280     rDY[1] = dY[1];
00281     rDY[2] = dY[2];
00282     rDY[3] = dY[3];
00283     rDY[4] = dY[4];
00284     rDY[5] = dY[5];
00285     rDY[6] = dY[6];
00286     rDY[7] = dY[7];
00287     rDY[8] = dY[8];
00288     rDY[9] = dY[9];
00289     rDY[10]= dY[10];
00290     rDY[11] = dY[11];
00291     rDY[12] = dY[12];
00292     rDY[13] = dY[13];
00293     rDY[14] = dY[14];
00294     rDY[15] = dY[15];
00295     rDY[16] = dY[16];
00296     rDY[17] = dY[17];
00297     rDY[18] = dY[18];
00298 }
00299 
00300 
00301 double TenTusscher2006OdeSystem::GetIIonic()
00302 {
00303       //Vector for state variables//
00304       double Y[19];
00305 
00306     Y[0] =  mStateVariables[0];// L_t Ype_Ca_current_d_gate_d (dimensionless)
00307     Y[1] =  mStateVariables[1];// L_t Ype_Ca_current_f2_gate_f2 (dimensionless)
00308     Y[2] =  mStateVariables[2];// L_t Ype_Ca_current_fCass_gate_fCass (dimensionless)
00309     Y[3] =  mStateVariables[3];// L_t Ype_Ca_current_f_gate_f (dimensionless)
00310     Y[4] =  mStateVariables[4];// calcium_d Ynamics_Ca_SR (millimolar)
00311     Y[5] =  mStateVariables[5];// calcium_d Ynamics_Ca_i (millimolar)
00312     Y[6] =  mStateVariables[6];// calcium_d Ynamics_Ca_ss (millimolar)
00313     Y[7] =  mStateVariables[7];// calcium_d Ynamics_R_prime (dimensionless)
00314     Y[8] =  mStateVariables[8];// fast_sodium_current_h_gate_h (dimensionless)
00315     Y[9] =  mStateVariables[9];// fast_sodium_current_j_gate_j (dimensionless)
00316     Y[10] =  mStateVariables[10];// fast_sodium_current_m_gate_m (dimensionless)
00317     Y[11] =  mStateVariables[11];// membrane_V (millivolt)
00318     Y[12] =  mStateVariables[12];// potassium_d Ynamics_K_i (millimolar)
00319     Y[13] =  mStateVariables[13];// rapid_time_dependent_potassium_current_Xr1_gate_Xr1 (dimensionless)
00320     Y[14] =  mStateVariables[14];// rapid_time_dependent_potassium_current_Xr2_gate_Xr2 (dimensionless)
00321     Y[15] =  mStateVariables[15];// slow_time_dependent_potassium_current_Xs_gate_Xs (dimensionless)
00322     Y[16] =  mStateVariables[16];// sodium_d Ynamics_Na_i (millimolar)
00323     Y[17] =  mStateVariables[17];// transient_outward_current_r_gate_r (dimensionless)
00324     Y[18] =  mStateVariables[18];// transient_outward_current_s_gate_s (dimensionless)
00325 
00326    double L_type_Ca_current_i_CaL = L_type_Ca_current_g_CaL*Y[0]*Y[3]*Y[1]*Y[2]*4.0*(Y[11]-15.0)*pow(membrane_F, 2.0)/(membrane_R*membrane_T)*(0.25*Y[6]*exp(2.0*(Y[11]-15.0)*membrane_F/(membrane_R*membrane_T))-calcium_dynamics_Ca_o)/(exp(2.0*(Y[11]-15.0)*membrane_F/(membrane_R*membrane_T))-1.0);
00327    double reversal_potentials_E_Ca = 0.5*membrane_R*membrane_T/membrane_F*log(calcium_dynamics_Ca_o/Y[5]);
00328    double calcium_background_current_i_b_Ca = calcium_background_current_g_bca*(Y[11]-reversal_potentials_E_Ca);
00329    double calcium_pump_current_i_p_Ca = calcium_pump_current_g_pCa*Y[5]/(Y[5]+calcium_pump_current_K_pCa);
00330    double sodium_calcium_exchanger_current_i_NaCa = sodium_calcium_exchanger_current_K_NaCa*(exp(sodium_calcium_exchanger_current_gamma*Y[11]*membrane_F/(membrane_R*membrane_T))*pow(Y[16], 3.0)*calcium_dynamics_Ca_o-exp((sodium_calcium_exchanger_current_gamma-1.0)*Y[11]*membrane_F/(membrane_R*membrane_T))*pow(sodium_dynamics_Na_o, 3.0)*Y[5]*sodium_calcium_exchanger_current_alpha)/((pow(sodium_calcium_exchanger_current_Km_Nai, 3.0)+pow(sodium_dynamics_Na_o, 3.0))*(sodium_calcium_exchanger_current_Km_Ca+calcium_dynamics_Ca_o)*(1.0+sodium_calcium_exchanger_current_K_sat*exp((sodium_calcium_exchanger_current_gamma-1.0)*Y[11]*membrane_F/(membrane_R*membrane_T))));
00331    double reversal_potentials_E_Na = membrane_R*membrane_T/membrane_F*log(sodium_dynamics_Na_o/Y[16]);
00332    double fast_sodium_current_i_Na = fast_sodium_current_g_Na*pow(Y[10], 3.0)*Y[8]*Y[9]*(Y[11]-reversal_potentials_E_Na);
00333 
00334    double reversal_potentials_E_K = membrane_R*membrane_T/membrane_F*log(potassium_dynamics_K_o/Y[12]);
00335    double inward_rectifier_potassium_current_alpha_K1 = 0.1/(1.0+exp(0.06*(Y[11]-reversal_potentials_E_K-200.0)));
00336    double inward_rectifier_potassium_current_beta_K1 = (3.0*exp(0.0002*(Y[11]-reversal_potentials_E_K+100.0))+exp(0.1*(Y[11]-reversal_potentials_E_K-10.0)))/(1.0+exp(-0.5*(Y[11]-reversal_potentials_E_K)));
00337    double inward_rectifier_potassium_current_xK1_inf = inward_rectifier_potassium_current_alpha_K1/(inward_rectifier_potassium_current_alpha_K1+inward_rectifier_potassium_current_beta_K1);
00338    double inward_rectifier_potassium_current_i_K1 = inward_rectifier_potassium_current_g_K1*inward_rectifier_potassium_current_xK1_inf*(Y[11]-reversal_potentials_E_K);
00339 
00340    double transient_outward_current_i_to = mScaleFactorIto*transient_outward_current_g_to*Y[17]*Y[18]*(Y[11]-reversal_potentials_E_K);
00341    double rapid_time_dependent_potassium_current_i_Kr = mScaleFactorGkr*rapid_time_dependent_potassium_current_g_Kr*sqrt(potassium_dynamics_K_o/5.4)*Y[13]*Y[14]*(Y[11]-reversal_potentials_E_K);
00342    double reversal_potentials_E_Ks = membrane_R*membrane_T/membrane_F*log((potassium_dynamics_K_o+reversal_potentials_P_kna*sodium_dynamics_Na_o)/(Y[12]+reversal_potentials_P_kna*Y[16]));
00343    double slow_time_dependent_potassium_current_i_Ks = mScaleFactorGks*slow_time_dependent_potassium_current_g_Ks*pow(Y[15], 2.0)*(Y[11]-reversal_potentials_E_Ks);
00344    double sodium_potassium_pump_current_i_NaK = sodium_potassium_pump_current_P_NaK*potassium_dynamics_K_o/(potassium_dynamics_K_o+sodium_potassium_pump_current_K_mk)*Y[16]/(Y[16]+sodium_potassium_pump_current_K_mNa)/(1.0+0.1245*exp(-0.1*Y[11]*membrane_F/(membrane_R*membrane_T))+0.0353*exp(-Y[11]*membrane_F/(membrane_R*membrane_T)));
00345    double sodium_background_current_i_b_Na = sodium_background_current_g_bna*(Y[11]-reversal_potentials_E_Na);
00346    double potassium_pump_current_i_p_K = potassium_pump_current_g_pK*(Y[11]-reversal_potentials_E_K)/(1.0+exp((25.0-Y[11])/5.98));
00347 
00348     double i_ionic = inward_rectifier_potassium_current_i_K1+transient_outward_current_i_to+rapid_time_dependent_potassium_current_i_Kr+slow_time_dependent_potassium_current_i_Ks+L_type_Ca_current_i_CaL+sodium_potassium_pump_current_i_NaK+fast_sodium_current_i_Na+sodium_background_current_i_b_Na+sodium_calcium_exchanger_current_i_NaCa+calcium_background_current_i_b_Ca+potassium_pump_current_i_p_K+calcium_pump_current_i_p_Ca; /*this is in nA*/
00349 
00350     assert(!std::isnan(i_ionic));
00351 
00352     double i_ionic_in_microA_per_cm2=i_ionic*1.0;
00353     return i_ionic_in_microA_per_cm2;
00354 
00355     /*   i_ionic for this model is in pA/pF.
00356      *    Please note that in the mono/bidomain formulation, i_ionic needs to be in microA/cm2.
00357      *    We then need to divide by the cell capacitance.
00358      *    The cell capacitance of the tenTusscher model is
00359      *    2.0 uF/cm2 in the paper
00360      *    0.185 uF/cm2 in this code (membrane_C)
00361      *    1.0 uF/cm2 in the EvaluateRhsDerivatives method above
00362      *
00363      *    For consistency, we choose the last option.
00364      *    i_ion*pow(10,-6) will be in microA/pF.
00365      *    Cm*pow(10,6) will be in pF/cm2.
00366      *    i_ion*pow(10,-6)*Cm*pow(10,6) = i_ion*Cm is in microA/cm2, the correct units
00367      */
00368 }
00369 
00370 void TenTusscher2006OdeSystem::VerifyStateVariables()
00371 {
00372     const std::vector<double>& rY = rGetStateVariables();
00373 
00374       double Y[19];
00375 
00376     Y[0] = rY[0];//  L_t Ype_Ca_current_d_gate_d (dimensionless)
00377     Y[1] = rY[1];//  L_t Ype_Ca_current_f2_gate_f2 (dimensionless)
00378     Y[2] = rY[2];//  L_t Ype_Ca_current_fCass_gate_fCass (dimensionless)
00379     Y[3] = rY[3];//  L_t Ype_Ca_current_f_gate_f (dimensionless)
00380     Y[4] = rY[4];//  calcium_d Ynamics_Ca_SR (millimolar)
00381     Y[5] = rY[5];//  calcium_d Ynamics_Ca_i (millimolar)
00382     Y[6] = rY[6];//  calcium_d Ynamics_Ca_ss (millimolar)
00383     Y[7] = rY[7];//  calcium_d Ynamics_R_prime (dimensionless)
00384     Y[8] = rY[8];//  fast_sodium_current_h_gate_h (dimensionless)
00385     Y[9] = rY[9];//  fast_sodium_current_j_gate_j (dimensionless)
00386     Y[10] = rY[10];//fast_sodium_current_m_gate_m (dimensionless)
00387     Y[11] = rY[11];//membrane_V (millivolt)
00388     Y[12] = rY[12];//potassium_d Ynamics_K_i (millimolar)
00389     Y[13] = rY[13];//rapid_time_dependent_potassium_current_Xr1_gate_Xr1 (dimensionless)
00390     Y[14] = rY[14];//rapid_time_dependent_potassium_current_Xr2_gate_Xr2 (dimensionless)
00391     Y[15] = rY[15];//slow_time_dependent_potassium_current_Xs_gate_Xs (dimensionless)
00392     Y[16] = rY[16];//sodium_d Ynamics_Na_i (millimolar)
00393     Y[17] = rY[17];//transient_outward_current_r_gate_r (dimensionless)
00394     Y[18] = rY[18];//transient_outward_current_s_gate_s (dimensionless)
00395 
00396     #define COVERAGE_IGNORE
00397     if (!(0<=Y[0] && Y[0]<=1))
00398     {
00399         EXCEPTION(DumpState("d gate of L type calcium channel is out of range!"));
00400     }
00401     if (!(0<=Y[1] && Y[1]<=1))
00402     {
00403         EXCEPTION(DumpState("f2 gate of L type calcium channel is out of range!"));
00404     }
00405     if (!(0<=Y[2] && Y[2]<=1))
00406     {
00407         EXCEPTION(DumpState("fCa gate of L type calcium channel is out of range!"));
00408     }
00409     if (!(0<=Y[3] && Y[3]<=1))
00410     {
00411         EXCEPTION(DumpState("f gate of L type calcium channel is out of range!"));
00412     }
00413     if (!(0<=Y[4]))
00414     {
00415         EXCEPTION(DumpState("CaSR is negative!"));
00416     }
00417     if (!(0<=Y[5]))
00418     {
00419         EXCEPTION(DumpState("Cai is negative!"));
00420     }
00421     if (!(0<=Y[6]))
00422     {
00423         EXCEPTION(DumpState("CaSS is negative!"));
00424     }
00425     if (!(0<=Y[7] && Y[7]<=1))
00426     {
00427         EXCEPTION(DumpState("R gate is out of range!"));
00428     }
00429     if (!(0<=Y[8] && Y[8]<=1))
00430     {
00431         EXCEPTION(DumpState("h gate of Na channel is out of range!"));
00432     }
00433     if (!(0<=Y[9] && Y[9]<=1))
00434     {
00435         EXCEPTION(DumpState("j gate of Na channel is out of range!"));
00436     }
00437     if (!(0<=Y[10] && Y[10]<=1))
00438     {
00439         EXCEPTION(DumpState("m gate of Na channel is out of range!"));
00440     }
00441     if (!(-200<=Y[11] && Y[11]<=200))
00442     {
00443         EXCEPTION(DumpState("Vm is REALLY out of range!"));
00444     }
00445     if (!(0<=Y[12]))
00446     {
00447         EXCEPTION(DumpState("Ki is negative!"));
00448     }
00449     if (!(0<=Y[13] && Y[13]<=1))
00450     {
00451         EXCEPTION(DumpState("xr1 gate of IKR channel is out of range!"));
00452     }
00453     if (!(0<=Y[14] && Y[14]<=1))
00454     {
00455         EXCEPTION(DumpState("xr2 gate of IKR channel is out of range!"));
00456     }
00457     if (!(0<=Y[15] && Y[15]<=1))
00458     {
00459         EXCEPTION(DumpState("xs1 gate of IKS channel is out of range!"));
00460     }
00461     if (!(0<=Y[16]))
00462     {
00463         EXCEPTION(DumpState("Nai is negative!"));
00464     }
00465     if (!(0<=Y[17] && Y[17]<=1))
00466     {
00467         EXCEPTION(DumpState("r gate of Ito channel is out of range!"));
00468     }
00469     if (!(0<=Y[18] && Y[18]<=1))
00470     {
00471         EXCEPTION(DumpState("s gate of Ito channel is out of range!"));
00472     }
00473 
00474     #undef COVERAGE_IGNORE
00475 }
00476 
00477 template<>
00478 void OdeSystemInformation<TenTusscher2006OdeSystem>::Initialise(void)
00479 {
00480     //Initailisation values are for steady state pacing at 500 ms BCL, please note that values may differ for different BCLs
00481    this->mVariableNames.push_back("d_gate_L");
00482    this->mVariableUnits.push_back("");
00483    this->mInitialConditions.push_back(0.00003373);
00484 
00485    this->mVariableNames.push_back("f2_gate_L");
00486    this->mVariableUnits.push_back("");
00487    this->mInitialConditions.push_back(0.9755);
00488 
00489    this->mVariableNames.push_back("fca_gate_L");
00490    this->mVariableUnits.push_back("");
00491    this->mInitialConditions.push_back(0.9953);
00492 
00493    this->mVariableNames.push_back("f_gate_L");
00494    this->mVariableUnits.push_back("");
00495    this->mInitialConditions.push_back(0.7888);
00496 
00497    this->mVariableNames.push_back("CaSR");
00498    this->mVariableUnits.push_back("mM");
00499    this->mInitialConditions.push_back(3.64);
00500 
00501    this->mVariableNames.push_back("Cai");
00502    this->mVariableUnits.push_back("mM");
00503    this->mInitialConditions.push_back(0.000126);
00504 
00505    this->mVariableNames.push_back("CaSS");
00506    this->mVariableUnits.push_back("mM");
00507    this->mInitialConditions.push_back(0.00036);
00508 
00509    this->mVariableNames.push_back("R_gate");
00510    this->mVariableUnits.push_back("");
00511    this->mInitialConditions.push_back(0.9073);
00512 
00513    this->mVariableNames.push_back("h_gate_Na");
00514    this->mVariableUnits.push_back("");
00515    this->mInitialConditions.push_back(0.7444);
00516 
00517    this->mVariableNames.push_back("j_gate_Na");
00518    this->mVariableUnits.push_back("");
00519    this->mInitialConditions.push_back(0.7045);
00520 
00521    this->mVariableNames.push_back("m_gate_Na");
00522    this->mVariableUnits.push_back("");
00523    this->mInitialConditions.push_back(0.00172);
00524 
00525    this->mVariableNames.push_back("V");
00526    this->mVariableUnits.push_back("mV");
00527    this->mInitialConditions.push_back(-85.23);
00528 
00529    this->mVariableNames.push_back("Ki");
00530    this->mVariableUnits.push_back("mM");
00531    this->mInitialConditions.push_back(136.89);
00532 
00533    this->mVariableNames.push_back("xr1_gate");
00534    this->mVariableUnits.push_back("");
00535    this->mInitialConditions.push_back(0.00621);
00536 
00537    this->mVariableNames.push_back("xr2_gate");
00538    this->mVariableUnits.push_back("");
00539    this->mInitialConditions.push_back(0.4712);
00540 
00541    this->mVariableNames.push_back("xS_gate");
00542    this->mVariableUnits.push_back("");
00543    this->mInitialConditions.push_back(0.0095);
00544 
00545    this->mVariableNames.push_back("Nai");
00546    this->mVariableUnits.push_back("mM");
00547    this->mInitialConditions.push_back(8.604);
00548 
00549    this->mVariableNames.push_back("r_gate_Ito");
00550    this->mVariableUnits.push_back("");
00551    this->mInitialConditions.push_back(0.0000000242);
00552 
00553    this->mVariableNames.push_back("s_gate_Ito");
00554    this->mVariableUnits.push_back("");
00555    this->mInitialConditions.push_back(0.999998);
00556 
00557     this->mInitialised = true;
00558 }
00559 
00560 
00561 

Generated by  doxygen 1.6.2