CorriasBuistSMCModified.cpp

/*

Copyright (C) University of Oxford, 2005-2011

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 <cmath>
#include <cassert>
#include <memory>
#include "Exception.hpp"
#include "OdeSystemInformation.hpp"
#include "CorriasBuistSMCModified.hpp"
#include "HeartConfig.hpp"


    CorriasBuistSMCModified::CorriasBuistSMCModified(boost::shared_ptr<AbstractIvpOdeSolver> pSolver, boost::shared_ptr<AbstractStimulusFunction> pIntracellularStimulus)
        : AbstractCardiacCell(
                pSolver,
                14,
                0,
                pIntracellularStimulus)
    {
        mpSystemInfo = OdeSystemInformation<CorriasBuistSMCModified>::Instance();
        mScaleFactorCarbonMonoxide = 1.0; // initialise to 1 --> no effect
        mFakeIccStimulusPresent = true;//by default we use the fake ICC stimulus

        /* parameters */
        Cm = 77.0*1e-6;// 77 pF --> microF

        Asurf_in_cm_square = Cm / HeartConfig::Instance()->GetCapacitance();
        Asurf = Asurf_in_cm_square / 0.01;//cm2 --> mm2 /*cell surface area (mm2)*/

        VolCell   =      3.5e-6;      /*cell volume (mm3)*/
        hCa       =      2.014e-4;    /*conc for half inactivation of fCa */
        sCa       =      1.31e-05;    /*slope factor for inactivation of fCa */

        /* concentrations */
        Ki         =     164.0;       /*intra K conc (mM)*/
        Nai        =     10.0;        /*intra Na conc (mM)*/
        ACh        =     1e-05;       /*acetylcholine conc (mM)*/
        CaiRest    =     0.9e-04;     /*baseline Ca conc (mM)*/

        /* maximum conductances*/
        gLVA_max    =    2.33766E-05; /*max conductance of ILVA*/                 // (0.18 nS) * 1e-6 (mS/nS) / Asurf (mm2) = mS/mm2
        gCaL_max    =    0.008441558; /*max conductance of ICaL*/                 // (65.0 nS) * 1e-6 (mS/nS) / Asurf (mm2) = mS/mm2
        gBK_max     =    0.005935065; /*max conductance of IBK)*/                 // (45.7 nS) * 1e-6 (mS/nS) / Asurf (mm2) = mS/mm2
        gKb_max     =    1.87013E-06; /*max conductance of IKb*/                  // (0.0144 nS) * 1e-6 (mS/nS) / Asurf (mm2) = mS/mm2
        gKA_max     =    0.001168831; /*max conductance of IKA*/                  // (9.0  nS) * 1e-6 (mS/nS) / Asurf (mm2) = mS/mm2
        gKr_max     =    0.004545455; /*max conductance of IKr*/                  // (35.0 nS) * 1e-6 (mS/nS) / Asurf (mm2) = mS/mm2
        gNa_max     =    0.00038961;  /*max conductance of INa*/                  // (3.0  nS) * 1e-6 (mS/nS) / Asurf (mm2) = mS/mm2
        gnsCC_max   =    0.006493506; /*max conductance of InsCC*/                // (50.0 nS) * 1e-6 (mS/nS) / Asurf (mm2) = mS/mm2
        gcouple     =     1.3*1e-6 / Asurf; /* coupling conductance bewteen fake ICC and SMC*/        // 1.3 nS * 1e-6 (mS/nS) / Asurf (mm2) = mS/mm2

        JCaExt_max   =   0.31705;     /*max flux of CaSR (mM/ms)*/

        /* Temperature corrections */
        Q10Ca        =   2.1;         /*(dim)*/
        Q10K         =   1.5;         /*(dim)*/ //1.365
        Q10Na        =   2.45;        /*(dim)*/
        Texp         =   297.0;       /*(degK)*/

         Ca_o   =         2.5;         // mM
         K_o     =        7.0;         // mM
         Na_o    =        137.0;       // mM

        /* Nernst parameters */
         R        =       8314.4;      // pJ/nmol/K
         T        =       310.0;       // degK
         F        =       96484.6;     // nC/nmol
         FoRT     =       0.03743;     // 1/mV
         RToF     =       26.7137;     // mV

        T_correct_Ca = pow(Q10Ca,(T-Texp)/10.0);/*temperature correction for Ca (dim)*/
        T_correct_K = pow(Q10K,(T-Texp)/10.0);  /*temperature correction for K (dim)*/
        T_correct_Na = pow(Q10Na,(T-Texp)/10.0);/*temperature correction for Na (dim)*/
        T_correct_gBK = gBK_max + 1.1*(T-Texp)*1e-6/Asurf; /*temperature correction for gBK*/  // (nS) * 1e-6 (mS/nS) / Asurf (mm2) = mS/mm2

        /* Nernst potentials */
        EK = RToF*log(K_o/Ki);                  /*Nernst potential for K (mV)*/
        ENa = RToF*log(Na_o/Nai);               /*Nernst potential for Na (mV)*/
        EnsCC = -28.0;                          /*Nernst potential for nsCC (mV)*/

        Init();

    }

    CorriasBuistSMCModified::~CorriasBuistSMCModified()
    {
    }

    void CorriasBuistSMCModified::VerifyStateVariables()
    {}

    void CorriasBuistSMCModified::SetCarbonMonoxideScaleFactor(double scaleFactor)
    {
        mScaleFactorCarbonMonoxide = scaleFactor;
    }

    void CorriasBuistSMCModified::SetFakeIccStimulusPresent(bool present)
    {
        mFakeIccStimulusPresent = present;
    }

    bool CorriasBuistSMCModified::GetFakeIccStimulusPresent()
    {
        return mFakeIccStimulusPresent;
    }

    double  CorriasBuistSMCModified::GetCarbonMonoxideScaleFactor()
    {
        return mScaleFactorCarbonMonoxide;
    }

    double CorriasBuistSMCModified::GetIIonic(const std::vector<double>* pStateVariables)
    {
        if (!pStateVariables) pStateVariables = &rGetStateVariables();
        const std::vector<double>& rY = *pStateVariables;

        double ECa = 0.5*RToF*log(Ca_o/rY[13]);

        /* inward sodium current */
        double INa = gNa_max*rY[6]*rY[7]*(rY[0]-ENa);

        /* L-type calcium current */
        double ICaL = gCaL_max*rY[1]*rY[2]*rY[3]*(rY[0]-ECa);

        /* low voltage activated (T-type) calcium current */
        double ILVA = gLVA_max*rY[4]*rY[5]*(rY[0]-ECa);

        /* large conductance calcium activated potassium current */
        double Po_BK = 1.0/(1.0+exp(-(rY[0]/17.0)-2.0*log(rY[13]/0.001)));
        double IBK = T_correct_gBK*Po_BK*(rY[0]-EK);

        /* delayed rectifier potassium current */
        double IKr = mScaleFactorCarbonMonoxide*gKr_max*rY[9]*rY[10]*(rY[0]-EK);

        /* A-type potassium current */
        double IKA = mScaleFactorCarbonMonoxide*gKA_max*rY[11]*rY[12]*(rY[0]-EK);

        /* background (leakage) potassium current */
        double IKb = mScaleFactorCarbonMonoxide*gKb_max*(rY[0]-EK);

        /* non-specific cation current */
        double hCa_nsCC = 1.0/(1.0+pow((rY[13]/0.0002),-4.0));
        double rACh_nsCC = 1.0/(1.0+(0.01/ACh));
        double InsCC = gnsCC_max*rY[8]*rACh_nsCC*hCa_nsCC*(rY[0]-EnsCC);

        /* phenomenological calcium extrusion current */
        double JCaExt = JCaExt_max*pow(rY[13],1.34);

        //i_ionic_in microA/mm2
        double i_ionic = INa+ICaL+ILVA+IKr+IKA+IBK+IKb+InsCC+(JCaExt*2.0*F*VolCell/Asurf);

        assert(!std::isnan(i_ionic));
        return i_ionic / 0.01;
    }

    void CorriasBuistSMCModified::EvaluateYDerivatives(double time, const std::vector<double>& rY, std::vector<double>& rDY)
    {

        // index [0] = -69.8;           // Vm       (mV)
        // index [1] = 5.0e-6;          // d_CaL    (dim)
        // index [2] = 0.953;           // f_CaL    (dim)
        // index [3] = 1.0;             // fCa_CaL  (dim)
        // index [4] = 0.0202;          // d_LVA    (dim)
        // index [5] = 1.0;             // f_LVA    (dim)
        // index [6] = 0.0086;          // d_Na     (dim)
        // index [7] = 0.061;           // f_Na     (dim)
        // index [8] = 0.096;           // m_nsCC   (dim)
        // index [9] = 1.92e-4;         // xr1      (dim)
        // index [10] = 0.812;          // xr2      (dim)
        // index [11] = 0.00415;        // xa1      (dim)
        // index [12] = 0.716;          // xa2      (dim)
        // index [13] = 0.9e-04;        // Cai      (mM)

        double ECa = 0.5*RToF*log(Ca_o/rY[13]);

        /* inward sodium current */
        double inf_d_Na = 1.0/(1.0+exp(-(rY[0]+47.0)/4.8));
        double tau_d_Na = (0.44-0.017*rY[0])*T_correct_Na;
        double inf_f_Na = 1.0/(1.0+exp((rY[0]+78.0)/3.0));
        double tau_f_Na = (5.5-0.25*rY[0])*T_correct_Na;

        double INa = gNa_max*rY[6]*rY[7]*(rY[0]-ENa);

        /* L-type calcium current */
        double inf_d_CaL = 1.0/(1.0+exp(-(rY[0]+17.0)/4.3));
        double tau_d_CaL = 0.47*T_correct_Ca;

        double inf_f_CaL = 1.0/(1.0+exp((rY[0]+43.0)/8.9));
        double tau_f_CaL = 86.0*T_correct_Ca;

        double inf_fCa_CaL = 1.0-(1.0/(1.0+exp(-((rY[13]-CaiRest)-hCa)/sCa)));
        double tau_fCa_CaL = 2.0*T_correct_Ca;
        double ICaL = gCaL_max*rY[1]*rY[2]*rY[3]*(rY[0]-ECa);

        /* low voltage activated (T-type) calcium current */
        double inf_d_LVA = 1.0/(1.0+exp(-(rY[0]+27.5)/10.9));
        double tau_d_LVA = 3.0*T_correct_Ca;

        double inf_f_LVA = 1.0/(1.0+exp((rY[0]+15.8)/7.0));
        double tau_f_LVA = 7.58*exp(rY[0]*0.00817)*T_correct_Ca;

        double ILVA = gLVA_max*rY[4]*rY[5]*(rY[0]-ECa);

        /* large conductance calcium activated potassium current */
        double Po_BK = 1.0/(1.0+exp(-(rY[0]/17.0)-2.0*log(rY[13]/0.001)));
        double IBK = T_correct_gBK*Po_BK*(rY[0]-EK);

        /* delayed rectifier potassium current */
        double inf_xr1 = 1.0/(1.0+exp(-(rY[0]+27.0)/5.0));
        double tau_xr1 = 80.0*T_correct_K;

        double inf_xr2 = 0.2+0.8/(1.0+exp((rY[0]+58.0)/10.0));
        double tau_xr2 = (-707.0+1481.0*exp((rY[0]+36.0)/95.0))*T_correct_K;

        double IKr = mScaleFactorCarbonMonoxide*gKr_max*rY[9]*rY[10]*(rY[0]-EK);

        /* A-type potassium current */
        double inf_xa1 = 1.0/(1.0+exp(-(rY[0]+26.5)/7.9));
        double tau_xa1 = (31.8+175.0*exp(-0.5*pow(((rY[0]+44.4)/22.3),2.0)))*T_correct_K;

        double inf_xa2 = 0.1+0.9/(1.0+exp((rY[0]+65.0)/6.2));
        double tau_xa2 = 90.0*T_correct_K;

        double IKA = mScaleFactorCarbonMonoxide*gKA_max*rY[11]*rY[12]*(rY[0]-EK);

        /* background (leakage) potassium current */
        double IKb = mScaleFactorCarbonMonoxide*gKb_max*(rY[0]-EK);

        /* non-specific cation current */
        double inf_m_nsCC = 1.0/(1.0+exp(-(rY[0]+25.0)/20.0));
        double tau_m_nsCC = 150.0/(1.0+exp(-(rY[0]+66.0)/26.0));
        double hCa_nsCC = 1.0/(1.0+pow((rY[13]/0.0002),-4.0));
        double rACh_nsCC = 1.0/(1.0+(0.01/ACh));

        double InsCC = gnsCC_max*rY[8]*rACh_nsCC*hCa_nsCC*(rY[0]-EnsCC);

        /* phenomenological calcium extrusion current */
        double JCaExt = JCaExt_max*pow(rY[13],1.34);


         double t_ICCplateau = 7582.0; // time_units
         double V_decay = 37.25; // voltage_units
         double t_ICCpeak = 98.0; // time_units

         double period = 20000.0; // time_units
         double stim_start = ((time > (period * 1.0)) && (time <= (period * 2.0))) ? (period * 1.0) : ((time > (period * 2.0)) && (time <= (period * 3.0))) ? (period * 2.0) : ((time > (period * 3.0)) && (time <= (period * 4.0))) ? (period * 3.0) : ((time > (period * 4.0)) && (time <= (period * 5.0))) ? (period * 4.0) : 0.0; // time_units
         double local_time = time - (stim_start + t_ICCpeak); // time_units
         double t_ICC_stimulus = 10000.0; // time_units
         double delta_VICC = 59.0; // voltage_units

        double i_stim;
        //see whether we are running this in isolaation (and we need the fake ICC stimulus) or coupled to a real ICC model
        if (mFakeIccStimulusPresent)
        {
            //for single cell simulations where we want the fake ICC stimulus in
            i_stim = (local_time < t_ICCpeak) ? (gcouple * delta_VICC) : ((local_time >= t_ICCpeak) && (local_time <= t_ICCplateau)) ? (gcouple * delta_VICC * (1.0 / (1.0 + exp((local_time - 8000.0) / 1000.0)))) : ((local_time > t_ICCplateau) && (local_time < t_ICC_stimulus)) ? (gcouple * V_decay * (1.0 / (1.0 + exp((local_time - 8000.0) / 150.0)))) : 0.0; // current_units
        }
        else
        {
            i_stim = GetStimulus(time);//for tissue simulations with current injected into SMC
        }

        /* membrane potential */
        double Iion = INa+ICaL+ILVA+IKr+IKA+IBK+IKb+InsCC+(JCaExt*2.0*F*VolCell/Asurf);

        double voltage_derivative;
        if (mSetVoltageDerivativeToZero)
        {
            voltage_derivative = 0.0;
        }
        else
        {
            voltage_derivative = (-1.0 / 0.01) * (-i_stim + Iion);//microA/mm2---> microA/cm2
            //std::cout<<rY[0]<<std::endl;
            assert(!std::isnan(voltage_derivative));
        }

        rDY[0] =  voltage_derivative;/* Vm */
        rDY[1] = (inf_d_CaL - rY[1])/tau_d_CaL;
        rDY[2] = (inf_f_CaL - rY[2])/tau_f_CaL;
        rDY[3] = (inf_fCa_CaL - rY[3])/tau_fCa_CaL;
        rDY[4] = (inf_d_LVA - rY[4])/tau_d_LVA;
        rDY[5] = (inf_f_LVA - rY[5])/tau_f_LVA;
        rDY[6] = (inf_d_Na - rY[6])/tau_d_Na;
        rDY[7] = (inf_f_Na - rY[7])/tau_f_Na;
        rDY[8] = (inf_m_nsCC - rY[8])/tau_m_nsCC;
        rDY[9] = (inf_xr1 - rY[9])/tau_xr1;
        rDY[10] = (inf_xr2 - rY[10])/tau_xr2;
        rDY[11] = (inf_xa1 - rY[11])/tau_xa1;
        rDY[12] = (inf_xa2 - rY[12])/tau_xa2;
        rDY[13] = (-(ICaL+ILVA)*Asurf/(2.0*F*VolCell)-JCaExt); /* intracellular calcium *1000 M-> mM; /1000 F units*/
    }

template<>
void OdeSystemInformation<CorriasBuistSMCModified>::Initialise(void)
{
    // Time units: time_units
    //
    this->mSystemName = "SMC_model_Martincode";

    this->mVariableNames.push_back("Vm_SM");
    this->mVariableUnits.push_back("mV");
    this->mInitialConditions.push_back(-69.8);           // Vm       (mV)

    this->mVariableNames.push_back("d_CaL");
    this->mVariableUnits.push_back("dim");
    this->mInitialConditions.push_back(5.0e-6);          // d_CaL    (dim));

    this->mVariableNames.push_back("f_CaL");
    this->mVariableUnits.push_back("dim");
    this->mInitialConditions.push_back(0.953);           // f_CaL    (dim)

    this->mVariableNames.push_back("fCa_CaL");
    this->mVariableUnits.push_back("dim");
    this->mInitialConditions.push_back(1.0);             // fCa_CaL  (dim)

    this->mVariableNames.push_back("d_LVA");
    this->mVariableUnits.push_back("dim");
    this->mInitialConditions.push_back(0.0202);          // d_LVA    (dim)

    this->mVariableNames.push_back("f_LVA");
    this->mVariableUnits.push_back("dim");
    this->mInitialConditions.push_back(1.0);             // f_LVA    (dim)

    this->mVariableNames.push_back("d_Na");
    this->mVariableUnits.push_back("dim");
    this->mInitialConditions.push_back(0.0086);          // d_Na     (dim)

    this->mVariableNames.push_back("f_Na");
    this->mVariableUnits.push_back("dim");
    this->mInitialConditions.push_back(0.061);           // f_Na     (dim)

    this->mVariableNames.push_back("m_nsCC");
    this->mVariableUnits.push_back("dim");
    this->mInitialConditions.push_back(0.096);           // m_nsCC   (dim)

    this->mVariableNames.push_back("xr1");
    this->mVariableUnits.push_back("dim");
    this->mInitialConditions.push_back(1.92e-4);         // xr1      (dim)

    this->mVariableNames.push_back("xr2");
    this->mVariableUnits.push_back("dim");
    this->mInitialConditions.push_back(0.812);          // xr2      (dim)

    this->mVariableNames.push_back("xa1");
    this->mVariableUnits.push_back("dim");
    this->mInitialConditions.push_back(0.00415);        // xa1      (dim)

    this->mVariableNames.push_back("xa2");
    this->mVariableUnits.push_back("dim");
    this->mInitialConditions.push_back(0.716);          // xa2      (dim)

    this->mVariableNames.push_back("Cai");
    this->mVariableUnits.push_back("mM");
    this->mInitialConditions.push_back(0.9e-04);        // Cai      (mM)

    this->mInitialised = true;
}


// Serialization for Boost >= 1.36
#include "SerializationExportWrapperForCpp.hpp"
CHASTE_CLASS_EXPORT(CorriasBuistSMCModified)