TenTusscher2006OdeSystem.cpp

/*

Copyright (C) University of Oxford, 2005-2009

University of Oxford means the Chancellor, Masters and Scholars of the
University of Oxford, having an administrative office at Wellington
Square, Oxford OX1 2JD, UK.

This file is part of Chaste.

Chaste is free software: you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License as published
by the Free Software Foundation, either version 2.1 of the License, or
(at your option) any later version.

Chaste is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
License for more details. The offer of Chaste under the terms of the
License is subject to the License being interpreted in accordance with
English Law and subject to any action against the University of Oxford
being under the jurisdiction of the English Courts.

You should have received a copy of the GNU Lesser General Public License
along with Chaste. If not, see <http://www.gnu.org/licenses/>.

*/
#include "TenTusscher2006OdeSystem.hpp"
#include "OdeSystemInformation.hpp"
#include <cmath>
#include <cstdio>
#include "Exception.hpp"


/*Constructor*/
TenTusscher2006OdeSystem::TenTusscher2006OdeSystem(
        AbstractIvpOdeSolver *pSolver,
        AbstractStimulusFunction *pIntracellularStimulus)
    : AbstractCardiacCell(pSolver, 19, 11, pIntracellularStimulus)
{

    mpSystemInfo = OdeSystemInformation<TenTusscher2006OdeSystem>::Instance();
    //Initialise the scale factors
    mScaleFactorGks=1.0;
    mScaleFactorIto=1.0;
    mScaleFactorGkr=1.0;
    
    AbstractCardiacCell::Init();
}

//Destructor
TenTusscher2006OdeSystem::~TenTusscher2006OdeSystem(void)
{
}

void TenTusscher2006OdeSystem::SetScaleFactorGks(double sfgks)
{
    mScaleFactorGks=sfgks;
}
void TenTusscher2006OdeSystem::SetScaleFactorIto(double sfito)
{
    mScaleFactorIto=sfito;
}
void TenTusscher2006OdeSystem::SetScaleFactorGkr(double sfgkr)
{
    mScaleFactorGkr=sfgkr;
}
void TenTusscher2006OdeSystem::EvaluateYDerivatives(double time,
                                                    const std::vector<double> &rY,
                                                    std::vector<double> &rDY)
{
   //---------------------------------------------------------------------------
   // State variables, initial value and names in comments beside
   //---------------------------------------------------------------------------
      //Vector for state variables//
      double Y[19];
      double dY[19];

    Y[0] = rY[0];// 3.373e-5;   // L_t Ype_Ca_current_d_gate_d (dimensionless)
    Y[1] = rY[1];// 0.9755;   // L_t Ype_Ca_current_f2_gate_f2 (dimensionless)
    Y[2] = rY[2];// 0.9953;   // L_t Ype_Ca_current_fCass_gate_fCass (dimensionless)
    Y[3] = rY[3];// 0.7888;   // L_t Ype_Ca_current_f_gate_f (dimensionless)
    Y[4] = rY[4];// 3.64;   // calcium_d Ynamics_Ca_SR (millimolar)
    Y[5] = rY[5];// 0.000126;   // calcium_d Ynamics_Ca_i (millimolar)
    Y[6] = rY[6];// 0.00036;   // calcium_d Ynamics_Ca_ss (millimolar)
    Y[7] = rY[7];// 0.9073;   // calcium_d Ynamics_R_prime (dimensionless)
    Y[8] = rY[8];// 0.7444;   // fast_sodium_current_h_gate_h (dimensionless)
    Y[9] = rY[9];// 0.7045;   // fast_sodium_current_j_gate_j (dimensionless)
    Y[10] = rY[10];// 0.00172;   // fast_sodium_current_m_gate_m (dimensionless)
    Y[11] = rY[11];// -85.23;   // membrane_V (millivolt)
    Y[12] = rY[12];// 136.89;   // potassium_d Ynamics_K_i (millimolar)
    Y[13] = rY[13];// 0.00621;   // rapid_time_dependent_potassium_current_Xr1_gate_Xr1 (dimensionless)
    Y[14] = rY[14];// 0.4712;   // rapid_time_dependent_potassium_current_Xr2_gate_Xr2 (dimensionless)
    Y[15] = rY[15];// 0.0095;   // slow_time_dependent_potassium_current_Xs_gate_Xs (dimensionless)
    Y[16] = rY[16];// 8.604;   // sodium_d Ynamics_Na_i (millimolar)
    Y[17] = rY[17];// 2.42e-8;   // transient_outward_current_r_gate_r (dimensionless)
    Y[18] = rY[18];// 0.999998;   // transient_outward_current_s_gate_s (dimensionless)

    VerifyStateVariables();

   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);
   L_type_Ca_current_d_gate_d_inf = 1.0/(1.0+exp((-8.0-Y[11])/7.5));
   L_type_Ca_current_d_gate_alpha_d = 1.4/(1.0+exp((-35.0-Y[11])/13.0))+0.25;
   L_type_Ca_current_d_gate_beta_d = 1.4/(1.0+exp((Y[11]+5.0)/5.0));
   L_type_Ca_current_d_gate_gamma_d = 1.0/(1.0+exp((50.0-Y[11])/20.0));
   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;
   dY[0] = (L_type_Ca_current_d_gate_d_inf-Y[0])/L_type_Ca_current_d_gate_tau_d;
   L_type_Ca_current_f2_gate_f2_inf = 0.67/(1.0+exp((Y[11]+35.0)/7.0))+0.33;
   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));
   dY[1] = (L_type_Ca_current_f2_gate_f2_inf-Y[1])/L_type_Ca_current_f2_gate_tau_f2;
   L_type_Ca_current_fCass_gate_fCass_inf = 0.6/(1.0+pow(Y[6]/0.05, 2.0))+0.4;
   L_type_Ca_current_fCass_gate_tau_fCass = 80.0/(1.0+pow(Y[6]/0.05, 2.0))+2.0;
   dY[2] = (L_type_Ca_current_fCass_gate_fCass_inf-Y[2])/L_type_Ca_current_fCass_gate_tau_fCass;
   L_type_Ca_current_f_gate_f_inf = 1.0/(1.0+exp((Y[11]+20.0)/7.0));
   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;
   dY[3] = (L_type_Ca_current_f_gate_f_inf-Y[3])/L_type_Ca_current_f_gate_tau_f;
   reversal_potentials_E_Ca = 0.5*membrane_R*membrane_T/membrane_F*log(calcium_dynamics_Ca_o/Y[5]);
   calcium_background_current_i_b_Ca = calcium_background_current_g_bca*(Y[11]-reversal_potentials_E_Ca);
   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));
   calcium_dynamics_k1 = calcium_dynamics_k1_prime/calcium_dynamics_kcasr;
   calcium_dynamics_O = calcium_dynamics_k1*pow(Y[6], 2.0)*Y[7]/(calcium_dynamics_k3+calcium_dynamics_k1*pow(Y[6], 2.0));
   calcium_dynamics_i_rel = calcium_dynamics_V_rel*calcium_dynamics_O*(Y[4]-Y[6]);
   calcium_dynamics_i_up = calcium_dynamics_Vmax_up/(1.0+pow(calcium_dynamics_K_up, 2.0)/pow(Y[5], 2.0));
   calcium_dynamics_i_leak = calcium_dynamics_V_leak*(Y[4]-Y[5]);
   calcium_dynamics_i_xfer = calcium_dynamics_V_xfer*(Y[6]-Y[5]);
   calcium_dynamics_k2 = calcium_dynamics_k2_prime*calcium_dynamics_kcasr;
   dY[7] = -calcium_dynamics_k2*Y[6]*Y[7]+calcium_dynamics_k4*(1.0-Y[7]);
   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));
   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));
   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));
   calcium_pump_current_i_p_Ca = calcium_pump_current_g_pCa*Y[5]/(Y[5]+calcium_pump_current_K_pCa);
   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))));
   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));
   dY[4] = calcium_dynamics_Ca_sr_bufsr*(calcium_dynamics_i_up-(calcium_dynamics_i_rel+calcium_dynamics_i_leak));
   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);
   reversal_potentials_E_Na = membrane_R*membrane_T/membrane_F*log(sodium_dynamics_Na_o/Y[16]);
   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);
   fast_sodium_current_h_gate_h_inf = 1.0/pow(1.0+exp((Y[11]+71.55)/7.43), 2.0);

   if (Y[11] < -40.0)
      fast_sodium_current_h_gate_alpha_h = 0.057*exp(-(Y[11]+80.0)/6.8);
   else
      fast_sodium_current_h_gate_alpha_h = 0.0;

   if (Y[11] < -40.0)
      fast_sodium_current_h_gate_beta_h = 2.7*exp(0.079*Y[11])+310000.0*exp(0.3485*Y[11]);
   else
      fast_sodium_current_h_gate_beta_h = 0.77/(0.13*(1.0+exp((Y[11]+10.66)/-11.1)));

   fast_sodium_current_h_gate_tau_h = 1.0/(fast_sodium_current_h_gate_alpha_h+fast_sodium_current_h_gate_beta_h);
   dY[8] = (fast_sodium_current_h_gate_h_inf-Y[8])/fast_sodium_current_h_gate_tau_h;
   fast_sodium_current_j_gate_j_inf = 1.0/pow(1.0+exp((Y[11]+71.55)/7.43), 2.0);

   if (Y[11] < -40.0)
      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)));
   else
      fast_sodium_current_j_gate_alpha_j = 0.0;

   if (Y[11] < -40.0)
      fast_sodium_current_j_gate_beta_j = 0.02424*exp(-0.01052*Y[11])/(1.0+exp(-0.1378*(Y[11]+40.14)));
   else
      fast_sodium_current_j_gate_beta_j = 0.6*exp(0.057*Y[11])/(1.0+exp(-0.1*(Y[11]+32.0)));

   fast_sodium_current_j_gate_tau_j = 1.0/(fast_sodium_current_j_gate_alpha_j+fast_sodium_current_j_gate_beta_j);
   dY[9] = (fast_sodium_current_j_gate_j_inf-Y[9])/fast_sodium_current_j_gate_tau_j;
   fast_sodium_current_m_gate_m_inf = 1.0/pow(1.0+exp((-56.86-Y[11])/9.03), 2.0);
   fast_sodium_current_m_gate_alpha_m = 1.0/(1.0+exp((-60.0-Y[11])/5.0));
   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));
   fast_sodium_current_m_gate_tau_m = fast_sodium_current_m_gate_alpha_m*fast_sodium_current_m_gate_beta_m;
   dY[10] = (fast_sodium_current_m_gate_m_inf-Y[10])/fast_sodium_current_m_gate_tau_m;
   reversal_potentials_E_K = membrane_R*membrane_T/membrane_F*log(potassium_dynamics_K_o/Y[12]);
   inward_rectifier_potassium_current_alpha_K1 = 0.1/(1.0+exp(0.06*(Y[11]-reversal_potentials_E_K-200.0)));
   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)));
   inward_rectifier_potassium_current_xK1_inf = inward_rectifier_potassium_current_alpha_K1/(inward_rectifier_potassium_current_alpha_K1+inward_rectifier_potassium_current_beta_K1);
   inward_rectifier_potassium_current_i_K1 = inward_rectifier_potassium_current_g_K1*inward_rectifier_potassium_current_xK1_inf*(Y[11]-reversal_potentials_E_K);


   transient_outward_current_i_to = mScaleFactorIto*transient_outward_current_g_to*Y[17]*Y[18]*(Y[11]-reversal_potentials_E_K);
   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);
   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]));
   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);
   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)));
   sodium_background_current_i_b_Na = sodium_background_current_g_bna*(Y[11]-reversal_potentials_E_Na);
   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));

   //stimulus current
    double i_stim = GetStimulus(time);

   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);
   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;
   rapid_time_dependent_potassium_current_Xr1_gate_xr1_inf = 1.0/(1.0+exp((-26.0-Y[11])/7.0));
   rapid_time_dependent_potassium_current_Xr1_gate_alpha_xr1 = 450.0/(1.0+exp((-45.0-Y[11])/10.0));
   rapid_time_dependent_potassium_current_Xr1_gate_beta_xr1 = 6.0/(1.0+exp((Y[11]+30.0)/11.5));
   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;
   dY[13] = (rapid_time_dependent_potassium_current_Xr1_gate_xr1_inf-Y[13])/rapid_time_dependent_potassium_current_Xr1_gate_tau_xr1;
   rapid_time_dependent_potassium_current_Xr2_gate_xr2_inf = 1.0/(1.0+exp((Y[11]+88.0)/24.0));
   rapid_time_dependent_potassium_current_Xr2_gate_alpha_xr2 = 3.0/(1.0+exp((-60.0-Y[11])/20.0));
   rapid_time_dependent_potassium_current_Xr2_gate_beta_xr2 = 1.12/(1.0+exp((Y[11]-60.0)/20.0));
   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;
   dY[14] = (rapid_time_dependent_potassium_current_Xr2_gate_xr2_inf-Y[14])/rapid_time_dependent_potassium_current_Xr2_gate_tau_xr2;
   slow_time_dependent_potassium_current_Xs_gate_xs_inf = 1.0/(1.0+exp((-5.0-Y[11])/14.0));
   slow_time_dependent_potassium_current_Xs_gate_alpha_xs = 1400.0/sqrt(1.0+exp((5.0-Y[11])/6.0));
   slow_time_dependent_potassium_current_Xs_gate_beta_xs = 1.0/(1.0+exp((Y[11]-35.0)/15.0));
   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;
   dY[15] = (slow_time_dependent_potassium_current_Xs_gate_xs_inf-Y[15])/slow_time_dependent_potassium_current_Xs_gate_tau_xs;
   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;
   transient_outward_current_r_gate_r_inf = 1.0/(1.0+exp((20.0-Y[11])/6.0));
   transient_outward_current_r_gate_tau_r = 9.5*exp(-pow(Y[11]+40.0, 2.0)/1800.0)+0.8;
   dY[17] = (transient_outward_current_r_gate_r_inf-Y[17])/transient_outward_current_r_gate_tau_r;
   transient_outward_current_s_gate_s_inf = 1.0/(1.0+exp((Y[11]+20.0)/5.0));
   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;
   dY[18] = (transient_outward_current_s_gate_s_inf-Y[18])/transient_outward_current_s_gate_tau_s;


    // do not update voltage if the mSetVoltageDerivativeToZero flag has been set
    if (mSetVoltageDerivativeToZero)
    {
        dY[11] = 0;
    }

    rDY[0] = dY[0];
    rDY[1] = dY[1];
    rDY[2] = dY[2];
    rDY[3] = dY[3];
    rDY[4] = dY[4];
    rDY[5] = dY[5];
    rDY[6] = dY[6];
    rDY[7] = dY[7];
    rDY[8] = dY[8];
    rDY[9] = dY[9];
    rDY[10]= dY[10];
    rDY[11] = dY[11];
    rDY[12] = dY[12];
    rDY[13] = dY[13];
    rDY[14] = dY[14];
    rDY[15] = dY[15];
    rDY[16] = dY[16];
    rDY[17] = dY[17];
    rDY[18] = dY[18];
}


double TenTusscher2006OdeSystem::GetIIonic()
{
      //Vector for state variables//
      double Y[19];
      double dY[19];
      
    Y[0] =  mStateVariables[0];// L_t Ype_Ca_current_d_gate_d (dimensionless)
    Y[1] =  mStateVariables[1];// L_t Ype_Ca_current_f2_gate_f2 (dimensionless)
    Y[2] =  mStateVariables[2];// L_t Ype_Ca_current_fCass_gate_fCass (dimensionless)
    Y[3] =  mStateVariables[3];// L_t Ype_Ca_current_f_gate_f (dimensionless)
    Y[4] =  mStateVariables[4];// calcium_d Ynamics_Ca_SR (millimolar)
    Y[5] =  mStateVariables[5];// calcium_d Ynamics_Ca_i (millimolar)
    Y[6] =  mStateVariables[6];// calcium_d Ynamics_Ca_ss (millimolar)
    Y[7] =  mStateVariables[7];// calcium_d Ynamics_R_prime (dimensionless)
    Y[8] =  mStateVariables[8];// fast_sodium_current_h_gate_h (dimensionless)
    Y[9] =  mStateVariables[9];// fast_sodium_current_j_gate_j (dimensionless)
    Y[10] =  mStateVariables[10];// fast_sodium_current_m_gate_m (dimensionless)
    Y[11] =  mStateVariables[11];// membrane_V (millivolt)
    Y[12] =  mStateVariables[12];// potassium_d Ynamics_K_i (millimolar)
    Y[13] =  mStateVariables[13];// rapid_time_dependent_potassium_current_Xr1_gate_Xr1 (dimensionless)
    Y[14] =  mStateVariables[14];// rapid_time_dependent_potassium_current_Xr2_gate_Xr2 (dimensionless)
    Y[15] =  mStateVariables[15];// slow_time_dependent_potassium_current_Xs_gate_Xs (dimensionless)
    Y[16] =  mStateVariables[16];// sodium_d Ynamics_Na_i (millimolar)
    Y[17] =  mStateVariables[17];// transient_outward_current_r_gate_r (dimensionless)
    Y[18] =  mStateVariables[18];// transient_outward_current_s_gate_s (dimensionless)
   
      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);
   L_type_Ca_current_d_gate_d_inf = 1.0/(1.0+exp((-8.0-Y[11])/7.5));
   L_type_Ca_current_d_gate_alpha_d = 1.4/(1.0+exp((-35.0-Y[11])/13.0))+0.25;
   L_type_Ca_current_d_gate_beta_d = 1.4/(1.0+exp((Y[11]+5.0)/5.0));
   L_type_Ca_current_d_gate_gamma_d = 1.0/(1.0+exp((50.0-Y[11])/20.0));
   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;
   dY[0] = (L_type_Ca_current_d_gate_d_inf-Y[0])/L_type_Ca_current_d_gate_tau_d;
   L_type_Ca_current_f2_gate_f2_inf = 0.67/(1.0+exp((Y[11]+35.0)/7.0))+0.33;
   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));
   dY[1] = (L_type_Ca_current_f2_gate_f2_inf-Y[1])/L_type_Ca_current_f2_gate_tau_f2;
   L_type_Ca_current_fCass_gate_fCass_inf = 0.6/(1.0+pow(Y[6]/0.05, 2.0))+0.4;
   L_type_Ca_current_fCass_gate_tau_fCass = 80.0/(1.0+pow(Y[6]/0.05, 2.0))+2.0;
   dY[2] = (L_type_Ca_current_fCass_gate_fCass_inf-Y[2])/L_type_Ca_current_fCass_gate_tau_fCass;
   L_type_Ca_current_f_gate_f_inf = 1.0/(1.0+exp((Y[11]+20.0)/7.0));
   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;
   dY[3] = (L_type_Ca_current_f_gate_f_inf-Y[3])/L_type_Ca_current_f_gate_tau_f;
   reversal_potentials_E_Ca = 0.5*membrane_R*membrane_T/membrane_F*log(calcium_dynamics_Ca_o/Y[5]);
   calcium_background_current_i_b_Ca = calcium_background_current_g_bca*(Y[11]-reversal_potentials_E_Ca);
   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));
   calcium_dynamics_k1 = calcium_dynamics_k1_prime/calcium_dynamics_kcasr;
   calcium_dynamics_O = calcium_dynamics_k1*pow(Y[6], 2.0)*Y[7]/(calcium_dynamics_k3+calcium_dynamics_k1*pow(Y[6], 2.0));
   calcium_dynamics_i_rel = calcium_dynamics_V_rel*calcium_dynamics_O*(Y[4]-Y[6]);
   calcium_dynamics_i_up = calcium_dynamics_Vmax_up/(1.0+pow(calcium_dynamics_K_up, 2.0)/pow(Y[5], 2.0));
   calcium_dynamics_i_leak = calcium_dynamics_V_leak*(Y[4]-Y[5]);
   calcium_dynamics_i_xfer = calcium_dynamics_V_xfer*(Y[6]-Y[5]);
   calcium_dynamics_k2 = calcium_dynamics_k2_prime*calcium_dynamics_kcasr;
   dY[7] = -calcium_dynamics_k2*Y[6]*Y[7]+calcium_dynamics_k4*(1.0-Y[7]);
   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));
   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));
   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));
   calcium_pump_current_i_p_Ca = calcium_pump_current_g_pCa*Y[5]/(Y[5]+calcium_pump_current_K_pCa);
   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))));
   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));
   dY[4] = calcium_dynamics_Ca_sr_bufsr*(calcium_dynamics_i_up-(calcium_dynamics_i_rel+calcium_dynamics_i_leak));
   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);
   reversal_potentials_E_Na = membrane_R*membrane_T/membrane_F*log(sodium_dynamics_Na_o/Y[16]);
   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);
   fast_sodium_current_h_gate_h_inf = 1.0/pow(1.0+exp((Y[11]+71.55)/7.43), 2.0);

   if (Y[11] < -40.0)
      fast_sodium_current_h_gate_alpha_h = 0.057*exp(-(Y[11]+80.0)/6.8);
   else
      fast_sodium_current_h_gate_alpha_h = 0.0;

   if (Y[11] < -40.0)
      fast_sodium_current_h_gate_beta_h = 2.7*exp(0.079*Y[11])+310000.0*exp(0.3485*Y[11]);
   else
      fast_sodium_current_h_gate_beta_h = 0.77/(0.13*(1.0+exp((Y[11]+10.66)/-11.1)));

   fast_sodium_current_h_gate_tau_h = 1.0/(fast_sodium_current_h_gate_alpha_h+fast_sodium_current_h_gate_beta_h);
   dY[8] = (fast_sodium_current_h_gate_h_inf-Y[8])/fast_sodium_current_h_gate_tau_h;
   fast_sodium_current_j_gate_j_inf = 1.0/pow(1.0+exp((Y[11]+71.55)/7.43), 2.0);

   if (Y[11] < -40.0)
      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)));
   else
      fast_sodium_current_j_gate_alpha_j = 0.0;

   if (Y[11] < -40.0)
      fast_sodium_current_j_gate_beta_j = 0.02424*exp(-0.01052*Y[11])/(1.0+exp(-0.1378*(Y[11]+40.14)));
   else
      fast_sodium_current_j_gate_beta_j = 0.6*exp(0.057*Y[11])/(1.0+exp(-0.1*(Y[11]+32.0)));

   fast_sodium_current_j_gate_tau_j = 1.0/(fast_sodium_current_j_gate_alpha_j+fast_sodium_current_j_gate_beta_j);
   dY[9] = (fast_sodium_current_j_gate_j_inf-Y[9])/fast_sodium_current_j_gate_tau_j;
   fast_sodium_current_m_gate_m_inf = 1.0/pow(1.0+exp((-56.86-Y[11])/9.03), 2.0);
   fast_sodium_current_m_gate_alpha_m = 1.0/(1.0+exp((-60.0-Y[11])/5.0));
   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));
   fast_sodium_current_m_gate_tau_m = fast_sodium_current_m_gate_alpha_m*fast_sodium_current_m_gate_beta_m;
   dY[10] = (fast_sodium_current_m_gate_m_inf-Y[10])/fast_sodium_current_m_gate_tau_m;
   reversal_potentials_E_K = membrane_R*membrane_T/membrane_F*log(potassium_dynamics_K_o/Y[12]);
   inward_rectifier_potassium_current_alpha_K1 = 0.1/(1.0+exp(0.06*(Y[11]-reversal_potentials_E_K-200.0)));
   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)));
   inward_rectifier_potassium_current_xK1_inf = inward_rectifier_potassium_current_alpha_K1/(inward_rectifier_potassium_current_alpha_K1+inward_rectifier_potassium_current_beta_K1);
   inward_rectifier_potassium_current_i_K1 = inward_rectifier_potassium_current_g_K1*inward_rectifier_potassium_current_xK1_inf*(Y[11]-reversal_potentials_E_K);


   transient_outward_current_i_to = mScaleFactorIto*transient_outward_current_g_to*Y[17]*Y[18]*(Y[11]-reversal_potentials_E_K);
   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);
   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]));
   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);
   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)));
   sodium_background_current_i_b_Na = sodium_background_current_g_bna*(Y[11]-reversal_potentials_E_Na);
   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));
    
   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);
   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;
   rapid_time_dependent_potassium_current_Xr1_gate_xr1_inf = 1.0/(1.0+exp((-26.0-Y[11])/7.0));
   rapid_time_dependent_potassium_current_Xr1_gate_alpha_xr1 = 450.0/(1.0+exp((-45.0-Y[11])/10.0));
   rapid_time_dependent_potassium_current_Xr1_gate_beta_xr1 = 6.0/(1.0+exp((Y[11]+30.0)/11.5));
   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;
   dY[13] = (rapid_time_dependent_potassium_current_Xr1_gate_xr1_inf-Y[13])/rapid_time_dependent_potassium_current_Xr1_gate_tau_xr1;
   rapid_time_dependent_potassium_current_Xr2_gate_xr2_inf = 1.0/(1.0+exp((Y[11]+88.0)/24.0));
   rapid_time_dependent_potassium_current_Xr2_gate_alpha_xr2 = 3.0/(1.0+exp((-60.0-Y[11])/20.0));
   rapid_time_dependent_potassium_current_Xr2_gate_beta_xr2 = 1.12/(1.0+exp((Y[11]-60.0)/20.0));
   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;
   dY[14] = (rapid_time_dependent_potassium_current_Xr2_gate_xr2_inf-Y[14])/rapid_time_dependent_potassium_current_Xr2_gate_tau_xr2;
   slow_time_dependent_potassium_current_Xs_gate_xs_inf = 1.0/(1.0+exp((-5.0-Y[11])/14.0));
   slow_time_dependent_potassium_current_Xs_gate_alpha_xs = 1400.0/sqrt(1.0+exp((5.0-Y[11])/6.0));
   slow_time_dependent_potassium_current_Xs_gate_beta_xs = 1.0/(1.0+exp((Y[11]-35.0)/15.0));
   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;
   dY[15] = (slow_time_dependent_potassium_current_Xs_gate_xs_inf-Y[15])/slow_time_dependent_potassium_current_Xs_gate_tau_xs;
   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;
   transient_outward_current_r_gate_r_inf = 1.0/(1.0+exp((20.0-Y[11])/6.0));
   transient_outward_current_r_gate_tau_r = 9.5*exp(-pow(Y[11]+40.0, 2.0)/1800.0)+0.8;
   dY[17] = (transient_outward_current_r_gate_r_inf-Y[17])/transient_outward_current_r_gate_tau_r;
   transient_outward_current_s_gate_s_inf = 1.0/(1.0+exp((Y[11]+20.0)/5.0));
   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;
   dY[18] = (transient_outward_current_s_gate_s_inf-Y[18])/transient_outward_current_s_gate_tau_s; 
    
    
    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*/
    
    assert(!isnan(i_ionic));

    double i_ionic_in_microA_per_cm2=i_ionic*membrane_Cm;
    return i_ionic_in_microA_per_cm2;
    
     /*   i_ionic for this model is in pA/pF. 
     *    Please note that in the mono/bidomain formulation, i_ionic needs to be in microA/cm2.
     *    The cell capacitance is, from the tenTusscher paper (code, actually), Cm = 0.185 microF/cm2.
     *    i_ion*pow(10,-6) will be in microA/pF. 
     *    Cm*pow(10,6) will be in pF/cm2.
     *    i_ion*pow(10,-6)*Cm*pow(10,6) = i_ion*Cm is in microA/cm2, the correct units
     */  

}

void TenTusscher2006OdeSystem::VerifyStateVariables()
{
    const std::vector<double>& rY = rGetStateVariables();

      double Y[19];

    Y[0] = rY[0];//  L_t Ype_Ca_current_d_gate_d (dimensionless)
    Y[1] = rY[1];//  L_t Ype_Ca_current_f2_gate_f2 (dimensionless)
    Y[2] = rY[2];//  L_t Ype_Ca_current_fCass_gate_fCass (dimensionless)
    Y[3] = rY[3];//  L_t Ype_Ca_current_f_gate_f (dimensionless)
    Y[4] = rY[4];//  calcium_d Ynamics_Ca_SR (millimolar)
    Y[5] = rY[5];//  calcium_d Ynamics_Ca_i (millimolar)
    Y[6] = rY[6];//  calcium_d Ynamics_Ca_ss (millimolar)
    Y[7] = rY[7];//  calcium_d Ynamics_R_prime (dimensionless)
    Y[8] = rY[8];//  fast_sodium_current_h_gate_h (dimensionless)
    Y[9] = rY[9];//  fast_sodium_current_j_gate_j (dimensionless)
    Y[10] = rY[10];//fast_sodium_current_m_gate_m (dimensionless)
    Y[11] = rY[11];//membrane_V (millivolt)
    Y[12] = rY[12];//potassium_d Ynamics_K_i (millimolar)
    Y[13] = rY[13];//rapid_time_dependent_potassium_current_Xr1_gate_Xr1 (dimensionless)
    Y[14] = rY[14];//rapid_time_dependent_potassium_current_Xr2_gate_Xr2 (dimensionless)
    Y[15] = rY[15];//slow_time_dependent_potassium_current_Xs_gate_Xs (dimensionless)
    Y[16] = rY[16];//sodium_d Ynamics_Na_i (millimolar)
    Y[17] = rY[17];//transient_outward_current_r_gate_r (dimensionless)
    Y[18] = rY[18];//transient_outward_current_s_gate_s (dimensionless)

    #define COVERAGE_IGNORE
    if (!(0<=Y[0] && Y[0]<=1))
    {
        EXCEPTION(DumpState("d gate of L type calcium channel is out of range!"));
    }
    if (!(0<=Y[1] && Y[1]<=1))
    {
        EXCEPTION(DumpState("f2 gate of L type calcium channel is out of range!"));
    }
    if (!(0<=Y[2] && Y[2]<=1))
    {
        EXCEPTION(DumpState("fCa gate of L type calcium channel is out of range!"));
    }
    if (!(0<=Y[3] && Y[3]<=1))
    {
        EXCEPTION(DumpState("f gate of L type calcium channel is out of range!"));
    }
    if (!(0<=Y[4]))
    {
        EXCEPTION(DumpState("CaSR is negative!"));
    }
    if (!(0<=Y[5]))
    {
        EXCEPTION(DumpState("Cai is negative!"));
    }
    if (!(0<=Y[6]))
    {
        EXCEPTION(DumpState("CaSS is negative!"));
    }
    if (!(0<=Y[7] && Y[7]<=1))
    {
        EXCEPTION(DumpState("R gate is out of range!"));
    }
    if (!(0<=Y[8] && Y[8]<=1))
    {
        EXCEPTION(DumpState("h gate of Na channel is out of range!"));
    }
    if (!(0<=Y[9] && Y[9]<=1))
    {
        EXCEPTION(DumpState("j gate of Na channel is out of range!"));
    }
    if (!(0<=Y[10] && Y[10]<=1))
    {
        EXCEPTION(DumpState("m gate of Na channel is out of range!"));
    }
    if (!(-200<=Y[11] && Y[11]<=200))
    {
        EXCEPTION(DumpState("Vm is REALLY out of range!"));
    }
    if (!(0<=Y[12]))
    {
        EXCEPTION(DumpState("Ki is negative!"));
    }
    if (!(0<=Y[13] && Y[13]<=1))
    {
        EXCEPTION(DumpState("xr1 gate of IKR channel is out of range!"));
    }
    if (!(0<=Y[14] && Y[14]<=1))
    {
        EXCEPTION(DumpState("xr2 gate of IKR channel is out of range!"));
    }
    if (!(0<=Y[15] && Y[15]<=1))
    {
        EXCEPTION(DumpState("xs1 gate of IKS channel is out of range!"));
    }
    if (!(0<=Y[16]))
    {
        EXCEPTION(DumpState("Nai is negative!"));
    }
    if (!(0<=Y[17] && Y[17]<=1))
    {
        EXCEPTION(DumpState("r gate of Ito channel is out of range!"));
    }
    if (!(0<=Y[18] && Y[18]<=1))
    {
        EXCEPTION(DumpState("s gate of Ito channel is out of range!"));
    }

    #undef COVERAGE_IGNORE
}

template<>
void OdeSystemInformation<TenTusscher2006OdeSystem>::Initialise(void)
{
    //Initailisation values are for steady state pacing at 500 ms BCL, please note that values may differ for different BCLs
   this->mVariableNames.push_back("d_gate_L");
   this->mVariableUnits.push_back("");
   this->mInitialConditions.push_back(0.00003373);

   this->mVariableNames.push_back("f2_gate_L");
   this->mVariableUnits.push_back("");
   this->mInitialConditions.push_back(0.9755);

   this->mVariableNames.push_back("fca_gate_L");
   this->mVariableUnits.push_back("");
   this->mInitialConditions.push_back(0.9953);

   this->mVariableNames.push_back("f_gate_L");
   this->mVariableUnits.push_back("");
   this->mInitialConditions.push_back(0.7888);

   this->mVariableNames.push_back("CaSR");
   this->mVariableUnits.push_back("mM");
   this->mInitialConditions.push_back(3.64);

   this->mVariableNames.push_back("Cai");
   this->mVariableUnits.push_back("mM");
   this->mInitialConditions.push_back(0.000126);

   this->mVariableNames.push_back("CaSS");
   this->mVariableUnits.push_back("mM");
   this->mInitialConditions.push_back(0.00036);

   this->mVariableNames.push_back("R_gate");
   this->mVariableUnits.push_back("");
   this->mInitialConditions.push_back(0.9073);

   this->mVariableNames.push_back("h_gate_Na");
   this->mVariableUnits.push_back("");
   this->mInitialConditions.push_back(0.7444);

   this->mVariableNames.push_back("j_gate_Na");
   this->mVariableUnits.push_back("");
   this->mInitialConditions.push_back(0.7045);

   this->mVariableNames.push_back("m_gate_Na");
   this->mVariableUnits.push_back("");
   this->mInitialConditions.push_back(0.00172);

   this->mVariableNames.push_back("V");
   this->mVariableUnits.push_back("mV");
   this->mInitialConditions.push_back(-85.23);

   this->mVariableNames.push_back("Ki");
   this->mVariableUnits.push_back("mM");
   this->mInitialConditions.push_back(136.89);

   this->mVariableNames.push_back("xr1_gate");
   this->mVariableUnits.push_back("");
   this->mInitialConditions.push_back(0.00621);

   this->mVariableNames.push_back("xr2_gate");
   this->mVariableUnits.push_back("");
   this->mInitialConditions.push_back(0.4712);

   this->mVariableNames.push_back("xS_gate");
   this->mVariableUnits.push_back("");
   this->mInitialConditions.push_back(0.0095);

   this->mVariableNames.push_back("Nai");
   this->mVariableUnits.push_back("mM");
   this->mInitialConditions.push_back(8.604);

   this->mVariableNames.push_back("r_gate_Ito");
   this->mVariableUnits.push_back("");
   this->mInitialConditions.push_back(0.0000000242);

   this->mVariableNames.push_back("s_gate_Ito");
   this->mVariableUnits.push_back("");
   this->mInitialConditions.push_back(0.999998);

    this->mInitialised = true;
}

---