CorriasBuistSMCModified.cpp

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