luo_rudy_1991CvodeOpt.cpp
00001 #ifdef CHASTE_CVODE
00013
00014 #include "luo_rudy_1991CvodeOpt.hpp"
00015 #include <cmath>
00016 #include <cassert>
00017 #include "Exception.hpp"
00018 #include "OdeSystemInformation.hpp"
00019
00020 class Cellluo_rudy_1991FromCellMLCvodeOpt_LookupTables
00021 {
00022 public:
00023 static Cellluo_rudy_1991FromCellMLCvodeOpt_LookupTables* Instance()
00024 {
00025 if (mpInstance == NULL)
00026 {
00027 mpInstance = new Cellluo_rudy_1991FromCellMLCvodeOpt_LookupTables;
00028 }
00029 return mpInstance;
00030 }
00031
00032
00033
00034 double* _lookup_0_row(unsigned i, double factor)
00035 {
00036 for (unsigned j=0; j<16; j++)
00037 {
00038 double y1 = _lookup_table_0[i][j];
00039 double y2 = _lookup_table_0[i+1][j];
00040 _lookup_table_0_row[j] = y1 + (y2-y1)*factor;
00041 }
00042 return _lookup_table_0_row;
00043 }
00044
00045
00046 protected:
00047 Cellluo_rudy_1991FromCellMLCvodeOpt_LookupTables(const Cellluo_rudy_1991FromCellMLCvodeOpt_LookupTables&);
00048 Cellluo_rudy_1991FromCellMLCvodeOpt_LookupTables& operator= (const Cellluo_rudy_1991FromCellMLCvodeOpt_LookupTables&);
00049 Cellluo_rudy_1991FromCellMLCvodeOpt_LookupTables()
00050 {
00051 assert(mpInstance == NULL);
00052 for (int i=0 ; i<16001; i++)
00053 {
00054 double var_membrane__V = -100.0001 + i*0.01;
00055 _lookup_table_0[i][0] = (0.32 * (var_membrane__V + 47.13)) / (1.0 - exp( -0.1 * (var_membrane__V + 47.13)));
00056 }
00057
00058 for (int i=0 ; i<16001; i++)
00059 {
00060 double var_membrane__V = -100.0001 + i*0.01;
00061 _lookup_table_0[i][1] = 0.08 * exp((-var_membrane__V) * 0.0909090909091);
00062 }
00063
00064 for (int i=0 ; i<16001; i++)
00065 {
00066 double var_membrane__V = -100.0001 + i*0.01;
00067 _lookup_table_0[i][2] = (var_membrane__V < -40.0) ? (0.135 * exp((80.0 + var_membrane__V) * -0.147058823529)) : 0.0;
00068 }
00069
00070 for (int i=0 ; i<16001; i++)
00071 {
00072 double var_membrane__V = -100.0001 + i*0.01;
00073 _lookup_table_0[i][3] = (var_membrane__V < -40.0) ? ((3.56 * exp(0.079 * var_membrane__V)) + (310000.0 * exp(0.35 * var_membrane__V))) : (1.0 / (0.13 * (1.0 + exp((var_membrane__V + 10.66) * -0.0900900900901))));
00074 }
00075
00076 for (int i=0 ; i<16001; i++)
00077 {
00078 double var_membrane__V = -100.0001 + i*0.01;
00079 _lookup_table_0[i][4] = (var_membrane__V < -40.0) ? (((( -127140.0 * exp(0.2444 * var_membrane__V)) - (3.474e-05 * exp( -0.04391 * var_membrane__V))) * (var_membrane__V + 37.78)) / (1.0 + exp(0.311 * (var_membrane__V + 79.23)))) : 0.0;
00080 }
00081
00082 for (int i=0 ; i<16001; i++)
00083 {
00084 double var_membrane__V = -100.0001 + i*0.01;
00085 _lookup_table_0[i][5] = (var_membrane__V < -40.0) ? ((0.1212 * exp( -0.01052 * var_membrane__V)) / (1.0 + exp( -0.1378 * (var_membrane__V + 40.14)))) : ((0.3 * exp( -2.535e-07 * var_membrane__V)) / (1.0 + exp( -0.1 * (var_membrane__V + 32.0))));
00086 }
00087
00088 for (int i=0 ; i<16001; i++)
00089 {
00090 double var_membrane__V = -100.0001 + i*0.01;
00091 _lookup_table_0[i][6] = (0.095 * exp( -0.01 * (var_membrane__V - 5.0))) / (1.0 + exp( -0.072 * (var_membrane__V - 5.0)));
00092 }
00093
00094 for (int i=0 ; i<16001; i++)
00095 {
00096 double var_membrane__V = -100.0001 + i*0.01;
00097 _lookup_table_0[i][7] = (0.07 * exp( -0.017 * (var_membrane__V + 44.0))) / (1.0 + exp(0.05 * (var_membrane__V + 44.0)));
00098 }
00099
00100 for (int i=0 ; i<16001; i++)
00101 {
00102 double var_membrane__V = -100.0001 + i*0.01;
00103 _lookup_table_0[i][8] = (0.012 * exp( -0.008 * (var_membrane__V + 28.0))) / (1.0 + exp(0.15 * (var_membrane__V + 28.0)));
00104 }
00105
00106 for (int i=0 ; i<16001; i++)
00107 {
00108 double var_membrane__V = -100.0001 + i*0.01;
00109 _lookup_table_0[i][9] = (0.0065 * exp( -0.02 * (var_membrane__V + 30.0))) / (1.0 + exp( -0.2 * (var_membrane__V + 30.0)));
00110 }
00111
00112 for (int i=0 ; i<16001; i++)
00113 {
00114 double var_membrane__V = -100.0001 + i*0.01;
00115 _lookup_table_0[i][10] = (var_membrane__V > -100.0) ? ((2.837 * (exp(0.04 * (var_membrane__V + 77.0)) - 1.0)) / ((var_membrane__V + 77.0) * exp(0.04 * (var_membrane__V + 35.0)))) : 1.0;
00116 }
00117
00118 for (int i=0 ; i<16001; i++)
00119 {
00120 double var_membrane__V = -100.0001 + i*0.01;
00121 _lookup_table_0[i][11] = (0.0005 * exp(0.083 * (var_membrane__V + 50.0))) / (1.0 + exp(0.057 * (var_membrane__V + 50.0)));
00122 }
00123
00124 for (int i=0 ; i<16001; i++)
00125 {
00126 double var_membrane__V = -100.0001 + i*0.01;
00127 _lookup_table_0[i][12] = (0.0013 * exp( -0.06 * (var_membrane__V + 20.0))) / (1.0 + exp( -0.04 * (var_membrane__V + 20.0)));
00128 }
00129
00130 for (int i=0 ; i<16001; i++)
00131 {
00132 double var_membrane__V = -100.0001 + i*0.01;
00133 _lookup_table_0[i][13] = ((0.49124 * exp(0.08032 * ((var_membrane__V + 5.476) - -87.8929017138))) + (1.0 * exp(0.06175 * (var_membrane__V - 506.417098286)))) / (1.0 + exp( -0.5143 * ((var_membrane__V - -87.8929017138) + 4.753)));
00134 }
00135
00136 for (int i=0 ; i<16001; i++)
00137 {
00138 double var_membrane__V = -100.0001 + i*0.01;
00139 _lookup_table_0[i][14] = 1.02 / (1.0 + exp(0.2385 * ((var_membrane__V - -87.8929017138) - 59.215)));
00140 }
00141
00142 for (int i=0 ; i<16001; i++)
00143 {
00144 double var_membrane__V = -100.0001 + i*0.01;
00145 _lookup_table_0[i][15] = 0.0183 * (1.0 / (1.0 + exp((7.488 - var_membrane__V) * 0.167224080268))) * (var_membrane__V - -87.8929017138);
00146 }
00147
00148 }
00149
00150 private:
00152 static Cellluo_rudy_1991FromCellMLCvodeOpt_LookupTables *mpInstance;
00153
00154
00155 double _lookup_table_0_row[16];
00156
00157
00158 double _lookup_table_0[16001][16];
00159
00160 };
00161
00162 Cellluo_rudy_1991FromCellMLCvodeOpt_LookupTables* Cellluo_rudy_1991FromCellMLCvodeOpt_LookupTables::mpInstance = NULL;
00163
00164 double Cellluo_rudy_1991FromCellMLCvodeOpt::Get_membrane__I_stim()
00165 {
00166 return var_membrane__I_stim;
00167 }
00168
00169 double Cellluo_rudy_1991FromCellMLCvodeOpt::Get_membrane__i_Na()
00170 {
00171 return var_membrane__i_Na;
00172 }
00173
00174 double Cellluo_rudy_1991FromCellMLCvodeOpt::Get_membrane__i_si()
00175 {
00176 return var_membrane__i_si;
00177 }
00178
00179 double Cellluo_rudy_1991FromCellMLCvodeOpt::Get_membrane__i_K()
00180 {
00181 return var_membrane__i_K;
00182 }
00183
00184 double Cellluo_rudy_1991FromCellMLCvodeOpt::Get_membrane__i_K1()
00185 {
00186 return var_membrane__i_K1;
00187 }
00188
00189 double Cellluo_rudy_1991FromCellMLCvodeOpt::Get_membrane__i_Kp()
00190 {
00191 return var_membrane__i_Kp;
00192 }
00193
00194 double Cellluo_rudy_1991FromCellMLCvodeOpt::Get_membrane__i_b()
00195 {
00196 return var_membrane__i_b;
00197 }
00198
00199 Cellluo_rudy_1991FromCellMLCvodeOpt::Cellluo_rudy_1991FromCellMLCvodeOpt(boost::shared_ptr<AbstractIvpOdeSolver> pOdeSolver , boost::shared_ptr<AbstractStimulusFunction> pIntracellularStimulus)
00200 : AbstractCvodeCell(
00201 pOdeSolver,
00202 8,
00203 0,
00204 pIntracellularStimulus)
00205 {
00206
00207
00208 mpSystemInfo = OdeSystemInformation<Cellluo_rudy_1991FromCellMLCvodeOpt>::Instance();
00209 Init();
00210
00211 }
00212
00213 Cellluo_rudy_1991FromCellMLCvodeOpt::~Cellluo_rudy_1991FromCellMLCvodeOpt()
00214 {
00215 }
00216
00217 void Cellluo_rudy_1991FromCellMLCvodeOpt::VerifyStateVariables()
00218 {}
00219
00220 double Cellluo_rudy_1991FromCellMLCvodeOpt::GetIIonic()
00221 {
00222 N_Vector rY = rGetStateVariables();
00223 double var_membrane__V = NV_Ith_S(rY, 0);
00224
00225 double var_fast_sodium_current_m_gate__m = NV_Ith_S(rY, 1);
00226
00227 double var_fast_sodium_current_h_gate__h = NV_Ith_S(rY, 2);
00228
00229 double var_fast_sodium_current_j_gate__j = NV_Ith_S(rY, 3);
00230
00231 double var_slow_inward_current_d_gate__d = NV_Ith_S(rY, 4);
00232
00233 double var_slow_inward_current_f_gate__f = NV_Ith_S(rY, 5);
00234
00235 double var_time_dependent_potassium_current_X_gate__X = NV_Ith_S(rY, 6);
00236
00237 double var_intracellular_calcium_concentration__Cai = NV_Ith_S(rY, 7);
00238
00239
00240
00241 if (var_membrane__V>59.9999 || var_membrane__V<-100.0001)
00242 {
00243 #define COVERAGE_IGNORE
00244 EXCEPTION(DumpState("V outside lookup table range", rY));
00245 #undef COVERAGE_IGNORE
00246 }
00247
00248 double _offset_0 = var_membrane__V - -100.0001;
00249 double _offset_0_over_table_step = _offset_0 * 100.0;
00250 unsigned _table_index_0 = (unsigned)(_offset_0_over_table_step);
00251 double _factor_0 = _offset_0_over_table_step - _table_index_0;
00252 double* _lt_0_row = Cellluo_rudy_1991FromCellMLCvodeOpt_LookupTables::Instance()->_lookup_0_row(_table_index_0, _factor_0);
00253
00254 var_membrane__i_Na = 23.0 * pow(var_fast_sodium_current_m_gate__m, 3.0) * var_fast_sodium_current_h_gate__h * var_fast_sodium_current_j_gate__j * (var_membrane__V - 54.7944639351);
00255 const double var_slow_inward_current__i_si = 0.09 * var_slow_inward_current_d_gate__d * var_slow_inward_current_f_gate__f * (var_membrane__V - (7.7 - (13.0287 * log(var_intracellular_calcium_concentration__Cai * 1.0))));
00256 var_membrane__i_si = var_slow_inward_current__i_si;
00257 var_membrane__i_K = 0.282 * var_time_dependent_potassium_current_X_gate__X * _lt_0_row[10] * (var_membrane__V - -77.5675843853);
00258 const double var_time_independent_potassium_current_K1_gate__alpha_K1 = _lt_0_row[14];
00259 var_membrane__i_K1 = 0.6047 * (var_time_independent_potassium_current_K1_gate__alpha_K1 / (var_time_independent_potassium_current_K1_gate__alpha_K1 + _lt_0_row[13])) * (var_membrane__V - -87.8929017138);
00260 var_membrane__i_Kp = _lt_0_row[15];
00261 var_membrane__i_b = 0.03921 * (var_membrane__V - -59.87);
00262
00263 return (var_membrane__i_Na+var_membrane__i_si+var_membrane__i_K+var_membrane__i_K1+var_membrane__i_Kp+var_membrane__i_b);
00264 }
00265
00266 void Cellluo_rudy_1991FromCellMLCvodeOpt::EvaluateRhs(double var_environment__time, const N_Vector rY, N_Vector rDY)
00267 {
00268
00269
00270 var_environment__time *= 1.0;
00271 double var_membrane__V = NV_Ith_S(rY, 0);
00272
00273 double var_fast_sodium_current_m_gate__m = NV_Ith_S(rY, 1);
00274
00275 double var_fast_sodium_current_h_gate__h = NV_Ith_S(rY, 2);
00276
00277 double var_fast_sodium_current_j_gate__j = NV_Ith_S(rY, 3);
00278
00279 double var_slow_inward_current_d_gate__d = NV_Ith_S(rY, 4);
00280
00281 double var_slow_inward_current_f_gate__f = NV_Ith_S(rY, 5);
00282
00283 double var_time_dependent_potassium_current_X_gate__X = NV_Ith_S(rY, 6);
00284
00285 double var_intracellular_calcium_concentration__Cai = NV_Ith_S(rY, 7);
00286
00287
00288
00289
00290 if (var_membrane__V>59.9999 || var_membrane__V<-100.0001)
00291 {
00292 #define COVERAGE_IGNORE
00293 EXCEPTION(DumpState("V outside lookup table range", rY));
00294 #undef COVERAGE_IGNORE
00295 }
00296
00297 double _offset_0 = var_membrane__V - -100.0001;
00298 double _offset_0_over_table_step = _offset_0 * 100.0;
00299 unsigned _table_index_0 = (unsigned)(_offset_0_over_table_step);
00300 double _factor_0 = _offset_0_over_table_step - _table_index_0;
00301 double* _lt_0_row = Cellluo_rudy_1991FromCellMLCvodeOpt_LookupTables::Instance()->_lookup_0_row(_table_index_0, _factor_0);
00302
00303
00304 var_membrane__I_stim = GetStimulus((1.0/1)*var_environment__time);
00305 var_membrane__i_Na = 23.0 * pow(var_fast_sodium_current_m_gate__m, 3.0) * var_fast_sodium_current_h_gate__h * var_fast_sodium_current_j_gate__j * (var_membrane__V - 54.7944639351);
00306 const double var_slow_inward_current__i_si = 0.09 * var_slow_inward_current_d_gate__d * var_slow_inward_current_f_gate__f * (var_membrane__V - (7.7 - (13.0287 * log(var_intracellular_calcium_concentration__Cai * 1.0))));
00307 var_membrane__i_si = var_slow_inward_current__i_si;
00308 var_membrane__i_K = 0.282 * var_time_dependent_potassium_current_X_gate__X * _lt_0_row[10] * (var_membrane__V - -77.5675843853);
00309 const double var_time_independent_potassium_current_K1_gate__alpha_K1 = _lt_0_row[14];
00310 var_membrane__i_K1 = 0.6047 * (var_time_independent_potassium_current_K1_gate__alpha_K1 / (var_time_independent_potassium_current_K1_gate__alpha_K1 + _lt_0_row[13])) * (var_membrane__V - -87.8929017138);
00311 var_membrane__i_Kp = _lt_0_row[15];
00312 var_membrane__i_b = 0.03921 * (var_membrane__V - -59.87);
00313
00314 double d_dt_membrane__V;
00315 if (mSetVoltageDerivativeToZero)
00316 {
00317 d_dt_membrane__V = 0.0;
00318 }
00319 else
00320 {
00321 d_dt_membrane__V = -1.0 * (var_membrane__I_stim + var_membrane__i_Na + var_membrane__i_si + var_membrane__i_K + var_membrane__i_K1 + var_membrane__i_Kp + var_membrane__i_b);
00322 }
00323
00324 const double d_dt_fast_sodium_current_m_gate__m = (_lt_0_row[0] * (1.0 - var_fast_sodium_current_m_gate__m)) - (_lt_0_row[1] * var_fast_sodium_current_m_gate__m);
00325 const double d_dt_fast_sodium_current_h_gate__h = (_lt_0_row[2] * (1.0 - var_fast_sodium_current_h_gate__h)) - (_lt_0_row[3] * var_fast_sodium_current_h_gate__h);
00326 const double d_dt_fast_sodium_current_j_gate__j = (_lt_0_row[4] * (1.0 - var_fast_sodium_current_j_gate__j)) - (_lt_0_row[5] * var_fast_sodium_current_j_gate__j);
00327 const double d_dt_slow_inward_current_d_gate__d = (_lt_0_row[6] * (1.0 - var_slow_inward_current_d_gate__d)) - (_lt_0_row[7] * var_slow_inward_current_d_gate__d);
00328 const double d_dt_slow_inward_current_f_gate__f = (_lt_0_row[8] * (1.0 - var_slow_inward_current_f_gate__f)) - (_lt_0_row[9] * var_slow_inward_current_f_gate__f);
00329 const double d_dt_time_dependent_potassium_current_X_gate__X = (_lt_0_row[11] * (1.0 - var_time_dependent_potassium_current_X_gate__X)) - (_lt_0_row[12] * var_time_dependent_potassium_current_X_gate__X);
00330 const double d_dt_intracellular_calcium_concentration__Cai = ( -0.0001 * var_slow_inward_current__i_si) + (0.07 * (0.0001 - var_intracellular_calcium_concentration__Cai));
00331
00332 NV_Ith_S(rDY, 0) = 1.0*d_dt_membrane__V;
00333 NV_Ith_S(rDY, 1) = 1.0*d_dt_fast_sodium_current_m_gate__m;
00334 NV_Ith_S(rDY, 2) = 1.0*d_dt_fast_sodium_current_h_gate__h;
00335 NV_Ith_S(rDY, 3) = 1.0*d_dt_fast_sodium_current_j_gate__j;
00336 NV_Ith_S(rDY, 4) = 1.0*d_dt_slow_inward_current_d_gate__d;
00337 NV_Ith_S(rDY, 5) = 1.0*d_dt_slow_inward_current_f_gate__f;
00338 NV_Ith_S(rDY, 6) = 1.0*d_dt_time_dependent_potassium_current_X_gate__X;
00339 NV_Ith_S(rDY, 7) = 1.0*d_dt_intracellular_calcium_concentration__Cai;
00340 }
00341
00342 template<>
00343 void OdeSystemInformation<Cellluo_rudy_1991FromCellMLCvodeOpt>::Initialise(void)
00344 {
00345
00346
00347 this->mVariableNames.push_back("membrane__V");
00348 this->mVariableUnits.push_back("millivolt");
00349 this->mInitialConditions.push_back(-83.853);
00350
00351 this->mVariableNames.push_back("fast_sodium_current_m_gate__m");
00352 this->mVariableUnits.push_back("dimensionless");
00353 this->mInitialConditions.push_back(0.00187018);
00354
00355 this->mVariableNames.push_back("fast_sodium_current_h_gate__h");
00356 this->mVariableUnits.push_back("dimensionless");
00357 this->mInitialConditions.push_back(0.9804713);
00358
00359 this->mVariableNames.push_back("fast_sodium_current_j_gate__j");
00360 this->mVariableUnits.push_back("dimensionless");
00361 this->mInitialConditions.push_back(0.98767124);
00362
00363 this->mVariableNames.push_back("slow_inward_current_d_gate__d");
00364 this->mVariableUnits.push_back("dimensionless");
00365 this->mInitialConditions.push_back(0.00316354);
00366
00367 this->mVariableNames.push_back("slow_inward_current_f_gate__f");
00368 this->mVariableUnits.push_back("dimensionless");
00369 this->mInitialConditions.push_back(0.99427859);
00370
00371 this->mVariableNames.push_back("time_dependent_potassium_current_X_gate__X");
00372 this->mVariableUnits.push_back("dimensionless");
00373 this->mInitialConditions.push_back(0.16647703);
00374
00375 this->mVariableNames.push_back("intracellular_calcium_concentration__Cai");
00376 this->mVariableUnits.push_back("millimolar");
00377 this->mInitialConditions.push_back(0.0002);
00378
00379 this->mInitialised = true;
00380 }
00381
00382
00383 #endif // CHASTE_CVODE
