CorriasBuistSMCModified.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
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 
00029 
00030 #include <cmath>
00031 #include <cassert>
00032 #include <memory>
00033 #include "Exception.hpp"
00034 #include "OdeSystemInformation.hpp"
00035 #include "CorriasBuistSMCModified.hpp"
00036 #include "HeartConfig.hpp"
00037 
00038 
00039     CorriasBuistSMCModified::CorriasBuistSMCModified(boost::shared_ptr<AbstractIvpOdeSolver> pSolver, boost::shared_ptr<AbstractStimulusFunction> pIntracellularStimulus)
00040         : AbstractCardiacCell(
00041                 pSolver,
00042                 14,
00043                 0,
00044                 pIntracellularStimulus)
00045     {
00046         mpSystemInfo = OdeSystemInformation<CorriasBuistSMCModified>::Instance();
00047         mScaleFactorCarbonMonoxide = 1.0; // initialise to 1 --> no effect
00048         mFakeIccStimulusPresent = true;//by default we use the fake ICC stimulus
00049 
00050         /* parameters */
00051         Cm = 77.0*1e-6;// 77 pF --> microF
00052 
00053         Asurf_in_cm_square = Cm / HeartConfig::Instance()->GetCapacitance();
00054         Asurf = Asurf_in_cm_square / 0.01;//cm2 --> mm2 /*cell surface area (mm2)*/
00055 
00056         VolCell   =      3.5e-6;      /*cell volume (mm3)*/
00057         hCa       =      2.014e-4;    /*conc for half inactivation of fCa */
00058         sCa       =      1.31e-05;    /*slope factor for inactivation of fCa */
00059 
00060         /* concentrations */
00061         Ki         =     164.0;       /*intra K conc (mM)*/
00062         Nai        =     10.0;        /*intra Na conc (mM)*/
00063         ACh        =     1e-05;       /*acetylcholine conc (mM)*/
00064         CaiRest    =     0.9e-04;     /*baseline Ca conc (mM)*/
00065 
00066         /* maximum conductances*/
00067         gLVA_max    =    2.33766E-05; /*max conductance of ILVA*/                 // (0.18 nS) * 1e-6 (mS/nS) / Asurf (mm2) = mS/mm2
00068         gCaL_max    =    0.008441558; /*max conductance of ICaL*/                 // (65.0 nS) * 1e-6 (mS/nS) / Asurf (mm2) = mS/mm2
00069         gBK_max     =    0.005935065; /*max conductance of IBK)*/                 // (45.7 nS) * 1e-6 (mS/nS) / Asurf (mm2) = mS/mm2
00070         gKb_max     =    1.87013E-06; /*max conductance of IKb*/                  // (0.0144 nS) * 1e-6 (mS/nS) / Asurf (mm2) = mS/mm2
00071         gKA_max     =    0.001168831; /*max conductance of IKA*/                  // (9.0  nS) * 1e-6 (mS/nS) / Asurf (mm2) = mS/mm2
00072         gKr_max     =    0.004545455; /*max conductance of IKr*/                  // (35.0 nS) * 1e-6 (mS/nS) / Asurf (mm2) = mS/mm2
00073         gNa_max     =    0.00038961;  /*max conductance of INa*/                  // (3.0  nS) * 1e-6 (mS/nS) / Asurf (mm2) = mS/mm2
00074         gnsCC_max   =    0.006493506; /*max conductance of InsCC*/                // (50.0 nS) * 1e-6 (mS/nS) / Asurf (mm2) = mS/mm2
00075         gcouple     =     1.3*1e-6 / Asurf; /* coupling conductance bewteen fake ICC and SMC*/        // 1.3 nS * 1e-6 (mS/nS) / Asurf (mm2) = mS/mm2
00076 
00077         JCaExt_max   =   0.31705;     /*max flux of CaSR (mM/ms)*/
00078 
00079         /* Temperature corrections */
00080         Q10Ca        =   2.1;         /*(dim)*/
00081         Q10K         =   1.5;         /*(dim)*/ //1.365
00082         Q10Na        =   2.45;        /*(dim)*/
00083         Texp         =   297.0;       /*(degK)*/
00084 
00085          Ca_o   =         2.5;         // mM
00086          K_o     =        7.0;         // mM
00087          Na_o    =        137.0;       // mM
00088 
00089         /* Nernst parameters */
00090          R        =       8314.4;      // pJ/nmol/K
00091          T        =       310.0;       // degK
00092          F        =       96484.6;     // nC/nmol
00093          FoRT     =       0.03743;     // 1/mV
00094          RToF     =       26.7137;     // mV
00095 
00096         T_correct_Ca = pow(Q10Ca,(T-Texp)/10.0);/*temperature correction for Ca (dim)*/
00097         T_correct_K = pow(Q10K,(T-Texp)/10.0);  /*temperature correction for K (dim)*/
00098         T_correct_Na = pow(Q10Na,(T-Texp)/10.0);/*temperature correction for Na (dim)*/
00099         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
00100 
00101         /* Nernst potentials */
00102         EK = RToF*log(K_o/Ki);                  /*Nernst potential for K (mV)*/
00103         ENa = RToF*log(Na_o/Nai);               /*Nernst potential for Na (mV)*/
00104         EnsCC = -28.0;                          /*Nernst potential for nsCC (mV)*/
00105 
00106         Init();
00107 
00108     }
00109 
00110     CorriasBuistSMCModified::~CorriasBuistSMCModified()
00111     {
00112     }
00113 
00114     void CorriasBuistSMCModified::VerifyStateVariables()
00115     {}
00116 
00117     void CorriasBuistSMCModified::SetCarbonMonoxideScaleFactor(double scaleFactor)
00118     {
00119         mScaleFactorCarbonMonoxide = scaleFactor;
00120     }
00121 
00122     void CorriasBuistSMCModified::SetFakeIccStimulusPresent(bool present)
00123     {
00124         mFakeIccStimulusPresent = present;
00125     }
00126 
00127     bool CorriasBuistSMCModified::GetFakeIccStimulusPresent()
00128     {
00129         return mFakeIccStimulusPresent;
00130     }
00131 
00132     double  CorriasBuistSMCModified::GetCarbonMonoxideScaleFactor()
00133     {
00134         return mScaleFactorCarbonMonoxide;
00135     }
00136 
00137     double CorriasBuistSMCModified::GetIIonic(const std::vector<double>* pStateVariables)
00138     {
00139         if (!pStateVariables) pStateVariables = &rGetStateVariables();
00140         const std::vector<double>& rY = *pStateVariables;
00141 
00142         double ECa = 0.5*RToF*log(Ca_o/rY[13]);
00143 
00144         /* inward sodium current */
00145         double INa = gNa_max*rY[6]*rY[7]*(rY[0]-ENa);
00146 
00147         /* L-type calcium current */
00148         double ICaL = gCaL_max*rY[1]*rY[2]*rY[3]*(rY[0]-ECa);
00149 
00150         /* low voltage activated (T-type) calcium current */
00151         double ILVA = gLVA_max*rY[4]*rY[5]*(rY[0]-ECa);
00152 
00153         /* large conductance calcium activated potassium current */
00154         double Po_BK = 1.0/(1.0+exp(-(rY[0]/17.0)-2.0*log(rY[13]/0.001)));
00155         double IBK = T_correct_gBK*Po_BK*(rY[0]-EK);
00156 
00157         /* delayed rectifier potassium current */
00158         double IKr = mScaleFactorCarbonMonoxide*gKr_max*rY[9]*rY[10]*(rY[0]-EK);
00159 
00160         /* A-type potassium current */
00161         double IKA = mScaleFactorCarbonMonoxide*gKA_max*rY[11]*rY[12]*(rY[0]-EK);
00162 
00163         /* background (leakage) potassium current */
00164         double IKb = mScaleFactorCarbonMonoxide*gKb_max*(rY[0]-EK);
00165 
00166         /* non-specific cation current */
00167         double hCa_nsCC = 1.0/(1.0+pow((rY[13]/0.0002),-4.0));
00168         double rACh_nsCC = 1.0/(1.0+(0.01/ACh));
00169         double InsCC = gnsCC_max*rY[8]*rACh_nsCC*hCa_nsCC*(rY[0]-EnsCC);
00170 
00171         /* phenomenological calcium extrusion current */
00172         double JCaExt = JCaExt_max*pow(rY[13],1.34);
00173 
00174         //i_ionic_in microA/mm2
00175         double i_ionic = INa+ICaL+ILVA+IKr+IKA+IBK+IKb+InsCC+(JCaExt*2.0*F*VolCell/Asurf);
00176 
00177         assert(!std::isnan(i_ionic));
00181         return i_ionic / 0.01;
00182     }
00183 
00184     void CorriasBuistSMCModified::EvaluateYDerivatives(double time, const std::vector<double>& rY, std::vector<double>& rDY)
00185     {
00186 
00187         // index [0] = -69.8;           // Vm       (mV)
00188         // index [1] = 5.0e-6;          // d_CaL    (dim)
00189         // index [2] = 0.953;           // f_CaL    (dim)
00190         // index [3] = 1.0;             // fCa_CaL  (dim)
00191         // index [4] = 0.0202;          // d_LVA    (dim)
00192         // index [5] = 1.0;             // f_LVA    (dim)
00193         // index [6] = 0.0086;          // d_Na     (dim)
00194         // index [7] = 0.061;           // f_Na     (dim)
00195         // index [8] = 0.096;           // m_nsCC   (dim)
00196         // index [9] = 1.92e-4;         // xr1      (dim)
00197         // index [10] = 0.812;          // xr2      (dim)
00198         // index [11] = 0.00415;        // xa1      (dim)
00199         // index [12] = 0.716;          // xa2      (dim)
00200         // index [13] = 0.9e-04;        // Cai      (mM)
00201 
00202         double ECa = 0.5*RToF*log(Ca_o/rY[13]);
00203 
00204         /* inward sodium current */
00205         double inf_d_Na = 1.0/(1.0+exp(-(rY[0]+47.0)/4.8));
00206         double tau_d_Na = (0.44-0.017*rY[0])*T_correct_Na;
00207         double inf_f_Na = 1.0/(1.0+exp((rY[0]+78.0)/3.0));
00208         double tau_f_Na = (5.5-0.25*rY[0])*T_correct_Na;
00209 
00210         double INa = gNa_max*rY[6]*rY[7]*(rY[0]-ENa);
00211 
00212         /* L-type calcium current */
00213         double inf_d_CaL = 1.0/(1.0+exp(-(rY[0]+17.0)/4.3));
00214         double tau_d_CaL = 0.47*T_correct_Ca;
00215 
00216         double inf_f_CaL = 1.0/(1.0+exp((rY[0]+43.0)/8.9));
00217         double tau_f_CaL = 86.0*T_correct_Ca;
00218 
00219         double inf_fCa_CaL = 1.0-(1.0/(1.0+exp(-((rY[13]-CaiRest)-hCa)/sCa)));
00220         double tau_fCa_CaL = 2.0*T_correct_Ca;
00221         double ICaL = gCaL_max*rY[1]*rY[2]*rY[3]*(rY[0]-ECa);
00222 
00223         /* low voltage activated (T-type) calcium current */
00224         double inf_d_LVA = 1.0/(1.0+exp(-(rY[0]+27.5)/10.9));
00225         double tau_d_LVA = 3.0*T_correct_Ca;
00226 
00227         double inf_f_LVA = 1.0/(1.0+exp((rY[0]+15.8)/7.0));
00228         double tau_f_LVA = 7.58*exp(rY[0]*0.00817)*T_correct_Ca;
00229 
00230         double ILVA = gLVA_max*rY[4]*rY[5]*(rY[0]-ECa);
00231 
00232         /* large conductance calcium activated potassium current */
00233         double Po_BK = 1.0/(1.0+exp(-(rY[0]/17.0)-2.0*log(rY[13]/0.001)));
00234         double IBK = T_correct_gBK*Po_BK*(rY[0]-EK);
00235 
00236         /* delayed rectifier potassium current */
00237         double inf_xr1 = 1.0/(1.0+exp(-(rY[0]+27.0)/5.0));
00238         double tau_xr1 = 80.0*T_correct_K;
00239 
00240         double inf_xr2 = 0.2+0.8/(1.0+exp((rY[0]+58.0)/10.0));
00241         double tau_xr2 = (-707.0+1481.0*exp((rY[0]+36.0)/95.0))*T_correct_K;
00242 
00243         double IKr = mScaleFactorCarbonMonoxide*gKr_max*rY[9]*rY[10]*(rY[0]-EK);
00244 
00245         /* A-type potassium current */
00246         double inf_xa1 = 1.0/(1.0+exp(-(rY[0]+26.5)/7.9));
00247         double tau_xa1 = (31.8+175.0*exp(-0.5*pow(((rY[0]+44.4)/22.3),2.0)))*T_correct_K;
00248 
00249         double inf_xa2 = 0.1+0.9/(1.0+exp((rY[0]+65.0)/6.2));
00250         double tau_xa2 = 90.0*T_correct_K;
00251 
00252         double IKA = mScaleFactorCarbonMonoxide*gKA_max*rY[11]*rY[12]*(rY[0]-EK);
00253 
00254         /* background (leakage) potassium current */
00255         double IKb = mScaleFactorCarbonMonoxide*gKb_max*(rY[0]-EK);
00256 
00257         /* non-specific cation current */
00258         double inf_m_nsCC = 1.0/(1.0+exp(-(rY[0]+25.0)/20.0));
00259         double tau_m_nsCC = 150.0/(1.0+exp(-(rY[0]+66.0)/26.0));
00260         double hCa_nsCC = 1.0/(1.0+pow((rY[13]/0.0002),-4.0));
00261         double rACh_nsCC = 1.0/(1.0+(0.01/ACh));
00262 
00263         double InsCC = gnsCC_max*rY[8]*rACh_nsCC*hCa_nsCC*(rY[0]-EnsCC);
00264 
00265         /* phenomenological calcium extrusion current */
00266         double JCaExt = JCaExt_max*pow(rY[13],1.34);
00267 
00268 
00269          double t_ICCplateau = 7582.0; // time_units
00270          double V_decay = 37.25; // voltage_units
00271          double t_ICCpeak = 98.0; // time_units
00272 
00273          double period = 20000.0; // time_units
00274          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
00275          double local_time = time - (stim_start + t_ICCpeak); // time_units
00276          double t_ICC_stimulus = 10000.0; // time_units
00277          double delta_VICC = 59.0; // voltage_units
00278 
00279         double i_stim;
00280         //see whether we are running this in isolaation (and we need the fake ICC stimulus) or coupled to a real ICC model
00281         if (mFakeIccStimulusPresent)
00282         {
00283             //for single cell simulations where we want the fake ICC stimulus in
00284             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
00285         }
00286         else
00287         {
00288             i_stim = GetStimulus(time);//for tissue simulations with current injected into SMC
00289         }
00290 
00291         /* membrane potential */
00292         double Iion = INa+ICaL+ILVA+IKr+IKA+IBK+IKb+InsCC+(JCaExt*2.0*F*VolCell/Asurf);
00293 
00294         double voltage_derivative;
00295         if (mSetVoltageDerivativeToZero)
00296         {
00297             voltage_derivative = 0.0;
00298         }
00299         else
00300         {
00301             voltage_derivative = (-1.0 / 0.01) * (-i_stim + Iion);//microA/mm2---> microA/cm2
00302             //std::cout<<rY[0]<<std::endl;
00303             assert(!std::isnan(voltage_derivative));
00304         }
00305 
00306         rDY[0] =  voltage_derivative;/* Vm */
00307         rDY[1] = (inf_d_CaL - rY[1])/tau_d_CaL;
00308         rDY[2] = (inf_f_CaL - rY[2])/tau_f_CaL;
00309         rDY[3] = (inf_fCa_CaL - rY[3])/tau_fCa_CaL;
00310         rDY[4] = (inf_d_LVA - rY[4])/tau_d_LVA;
00311         rDY[5] = (inf_f_LVA - rY[5])/tau_f_LVA;
00312         rDY[6] = (inf_d_Na - rY[6])/tau_d_Na;
00313         rDY[7] = (inf_f_Na - rY[7])/tau_f_Na;
00314         rDY[8] = (inf_m_nsCC - rY[8])/tau_m_nsCC;
00315         rDY[9] = (inf_xr1 - rY[9])/tau_xr1;
00316         rDY[10] = (inf_xr2 - rY[10])/tau_xr2;
00317         rDY[11] = (inf_xa1 - rY[11])/tau_xa1;
00318         rDY[12] = (inf_xa2 - rY[12])/tau_xa2;
00319         rDY[13] = (-(ICaL+ILVA)*Asurf/(2.0*F*VolCell)-JCaExt); /* intracellular calcium *1000 M-> mM; /1000 F units*/
00320     }
00321 
00322 template<>
00323 void OdeSystemInformation<CorriasBuistSMCModified>::Initialise(void)
00324 {
00325     // Time units: time_units
00326     //
00327     this->mSystemName = "SMC_model_Martincode";
00328 
00329     this->mVariableNames.push_back("Vm_SM");
00330     this->mVariableUnits.push_back("mV");
00331     this->mInitialConditions.push_back(-69.8);           // Vm       (mV)
00332 
00333     this->mVariableNames.push_back("d_CaL");
00334     this->mVariableUnits.push_back("dim");
00335     this->mInitialConditions.push_back(5.0e-6);          // d_CaL    (dim));
00336 
00337     this->mVariableNames.push_back("f_CaL");
00338     this->mVariableUnits.push_back("dim");
00339     this->mInitialConditions.push_back(0.953);           // f_CaL    (dim)
00340 
00341     this->mVariableNames.push_back("fCa_CaL");
00342     this->mVariableUnits.push_back("dim");
00343     this->mInitialConditions.push_back(1.0);             // fCa_CaL  (dim)
00344 
00345     this->mVariableNames.push_back("d_LVA");
00346     this->mVariableUnits.push_back("dim");
00347     this->mInitialConditions.push_back(0.0202);          // d_LVA    (dim)
00348 
00349     this->mVariableNames.push_back("f_LVA");
00350     this->mVariableUnits.push_back("dim");
00351     this->mInitialConditions.push_back(1.0);             // f_LVA    (dim)
00352 
00353     this->mVariableNames.push_back("d_Na");
00354     this->mVariableUnits.push_back("dim");
00355     this->mInitialConditions.push_back(0.0086);          // d_Na     (dim)
00356 
00357     this->mVariableNames.push_back("f_Na");
00358     this->mVariableUnits.push_back("dim");
00359     this->mInitialConditions.push_back(0.061);           // f_Na     (dim)
00360 
00361     this->mVariableNames.push_back("m_nsCC");
00362     this->mVariableUnits.push_back("dim");
00363     this->mInitialConditions.push_back(0.096);           // m_nsCC   (dim)
00364 
00365     this->mVariableNames.push_back("xr1");
00366     this->mVariableUnits.push_back("dim");
00367     this->mInitialConditions.push_back(1.92e-4);         // xr1      (dim)
00368 
00369     this->mVariableNames.push_back("xr2");
00370     this->mVariableUnits.push_back("dim");
00371     this->mInitialConditions.push_back(0.812);          // xr2      (dim)
00372 
00373     this->mVariableNames.push_back("xa1");
00374     this->mVariableUnits.push_back("dim");
00375     this->mInitialConditions.push_back(0.00415);        // xa1      (dim)
00376 
00377     this->mVariableNames.push_back("xa2");
00378     this->mVariableUnits.push_back("dim");
00379     this->mInitialConditions.push_back(0.716);          // xa2      (dim)
00380 
00381     this->mVariableNames.push_back("Cai");
00382     this->mVariableUnits.push_back("mM");
00383     this->mInitialConditions.push_back(0.9e-04);        // Cai      (mM)
00384 
00385     this->mInitialised = true;
00386 }
00387 
00388 
00389 // Serialization for Boost >= 1.36
00390 #include "SerializationExportWrapperForCpp.hpp"
00391 CHASTE_CLASS_EXPORT(CorriasBuistSMCModified)
Generated on Thu Dec 22 13:00:06 2011 for Chaste by  doxygen 1.6.3