- Generated on Mon Nov 5 2012 12:47:14 for Chaste by
1.7.4
|
Chaste Release::3.1
|
00001 /* 00002 00003 Copyright (c) 2005-2012, 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)
1.7.4