1
// Copyright (C) Glenn Terje Lines, Ola Skavhaug and Simula Research Laboratory.
2
// Licensed under the GNU LGPL Version 2.1.
4
// Original code copied from PyCC.
5
// Modified by Anders Logg 2006.
7
// First added: 2006-05-24
8
// Last changed: 2006-08-21
10
// This demo solves the Courtemanche model for cardiac excitation.
14
using namespace dolfin;
16
class Courtemanche : public ODE
20
Courtemanche() : ODE(21, 300.0)
45
I_NaKmax = 0.59933874;
77
message("Function evaluations: %d", num_fevals);
78
message("Potential at end time: %.6f", VT);
81
void u0(uBlasVector& u)
110
void f(const uBlasVector& u, real t, uBlasVector& y)
113
computeGateCoefficients(u);
115
y(0) = -1.0/Cm*(I_ion + ist);
116
y(1) = (m_inf - m)/tau_m;
117
y(2) = (h_inf - h)/tau_h;
118
y(3) = (j_inf - j)/tau_j;
119
y(4) = (oa_inf - oa)/tau_oa;
120
y(5) = (oi_inf - oi)/tau_oi;
121
y(6) = (ua_inf - ua)/tau_ua;
122
y(7) = (ui_inf - ui)/tau_ui;
123
y(8) = (xr_inf - xr)/tau_xr;
124
y(9) = (xs_inf - xs)/tau_xs;
125
y(10) = (d_inf - d)/tau_d;
126
y(11) = (f_inf - ff)/tau_f;
127
y(12) = (fca_inf - fca)/tau_fca;
128
y(13) = (u_inf - uu)/tau_u;
129
y(14) = (v_inf - v)/tau_v;
130
y(15) = (w_inf - w)/tau_w;
131
y(16) = (-3.0*I_NaK - 3.0*I_NaCa - I_bNa - I_Na)/(F*Vi);
133
y(18) = (I_tr - I_rel)/(1.0 + Csqn_max*K_mCsqn/((Ca_rel + K_mCsqn)*(Ca_rel + K_mCsqn)));
134
y(19) = I_up - I_upleak - I_tr*(Vrel/Vup);
135
y(20) = (2.0*I_NaK - I_K1 - I_to - I_Kur - I_Kr - I_Ks)/(F*Vi);
140
void computeCurrents(const uBlasVector& u)
164
I_rel = k_rel*uu*uu*v*w*(Ca_rel - Ca_i);
165
I_CaL = Cm*g_CaL*d*ff*fca*(V - 65);
166
I_NaCa = Cm*(I_NaCamax*(exp(gamma*F*V/(R*T))*Na_i*Na_i*Na_i*Ca_o - exp((gamma - 1)*F*V/(R*T))*Na_o*Na_o*Na_o*Ca_i))/((K_mNa*K_mNa*K_mNa + Na_o*Na_o*Na_o)*(K_mCa + Ca_o)*(1 + k_sat*exp((gamma -1 )*F*V/(R*T))));
167
//Fn = 1e-12*Vrel*I_rel - 5e-13/F*(0.5*I_CaL - 0.2*I_NaCa);
168
sigma = (1.0/7.0)*(exp(Na_o/67.3) - 1.0);
169
f_NaK = 1.0/(1.0 + 0.1245*exp(-0.1*F*V/(R*T))+ 0.0365*sigma*exp(-F*V/(R*T)));
170
I_NaK = Cm*I_NaKmax*f_NaK/(1.0 + pow((K_mNai/Na_i),1.5))*(K_o/(K_o + K_mKo));
171
E_Na = R*T/(z_Na*F)*log(Na_o/Na_i);
172
I_bNa = Cm*g_bNa*(V - E_Na);
173
I_Na = Cm*g_Na*m*m*m*h*j*(V - E_Na);
174
I_pCa = Cm*I_pCamax*Ca_i/(0.0005 + Ca_i);
175
E_Ca = R*T/(z_Ca*F)*log(Ca_o/Ca_i);
176
I_bCa = Cm*g_bCa*(V - E_Ca);
177
I_upleak = (Ca_up/Ca_upmax)*I_upmax;
178
I_up = I_upmax/(1.0 + (K_up/Ca_i));
179
I_tr = (Ca_up - Ca_rel)/tau_tr;
180
E_K = R*T/(z_K*F)*log(K_o/K_i);
181
I_K1 = Cm*(g_K1*(V - E_K))/(1.0 + exp(0.07*(V + 80.0)));
182
I_to = Cm*g_to*oa*oa*oa*oi*(V - E_K);
183
g_Kur = 0.005 + 0.05/(1.0 + exp((V - 15.0)/-13.0));
184
I_Kur = Cm*g_Kur*ua*ua*ua*ui*(V - E_K);
185
I_Kr = Cm*(g_Kr*xr*(V - E_K))/(1.0 + exp((V + 15.0)/22.4));
186
I_Ks = Cm*g_Ks*xs*xs*(V - E_K);
187
I_ion = I_Na + I_K1 + I_to + I_Kur + I_Kr + I_Ks + I_CaL + I_pCa + I_NaK + I_NaCa + I_bNa + I_bCa;
188
B1 = (2.0*I_NaCa - I_pCa - I_CaL - I_bCa)/(2.0*F*Vi) + (Vup*(I_upleak - I_up) + I_rel*Vrel)/Vi;
189
B2 = 1.0 + Trpn_max*K_mTrpn/((Ca_i + K_mTrpn)*(Ca_i + K_mTrpn)) + Cmdn_max*K_mCmdn/((Ca_i + K_mCmdn)*(Ca_i + K_mCmdn));
192
void computeGateCoefficients(const uBlasVector& u)
199
alpha_m = 0.32*(V + 47.13)/(1.0 - exp(-0.1*(V + 47.13)));
201
beta_m = 0.08*exp(V/-11.0);
202
tau_m = 1.0/(alpha_m + beta_m);
203
m_inf = alpha_m*tau_m;
206
beta_h = 1.0/(0.13*(1.0 + exp((V + 10.66)/-11.1)));
208
alpha_h = 0.135*exp((V + 80.0)/-6.8);
209
beta_h = 3.56*exp(0.079*V)+3.1e5*exp(0.35*V);
212
tau_h = 1.0/(alpha_h + beta_h);
213
h_inf = alpha_h*tau_h;
218
beta_j = 0.3*(exp(-2.535e-7*V))/(1.0 + exp(-0.1*(V +32.0)));
220
alpha_j = (-127140.0*exp(0.2444*V)-3.474e-5*exp(-0.04391*V))*(V + 37.78)/(1.0 + exp(0.311*(V + 79.23)));
221
beta_j = 0.1212*(exp(-0.01052*V))/(1.0 + exp(-0.1378*(V + 40.14)));
224
tau_j = 1.0/(alpha_j + beta_j);
225
j_inf = alpha_j*tau_j;
227
alpha_oa = 0.65/(exp((V + 10.0)/-8.5) + exp((V - 30.0)/-59.0));
228
beta_oa = 0.65/(2.5 + exp((V + 82.0)/17.0));
229
tau_oa = 1.0/((alpha_oa + beta_oa)*K_Q10);
230
oa_inf = 1.0/(1.0 + exp((V + 20.47)/-17.54));
232
alpha_oi = 1.0/(18.53 + exp((V + 113.7)/10.95));
233
beta_oi = 1.0/(35.56 + exp((V + 1.26)/-7.44));
234
tau_oi = 1.0/((alpha_oi + beta_oi)*K_Q10);
235
oi_inf = 1.0/(1.0 + exp((V + 43.1)/5.3));
237
alpha_ua = 0.65/(exp((V + 10.0)/-8.5) + exp((V - 30)/-59.0));
238
beta_ua = 0.65/(2.5 + exp((V + 82.0)/17.0));
239
tau_ua = 1.0/((alpha_ua + beta_ua)*K_Q10);
240
ua_inf = 1.0/(1.0 + exp((V + 30.3)/-9.6));
242
alpha_ui = 1.0/(21.0 + exp((V - 185.0)/-28.0));
243
beta_ui = exp((V - 158.0)/16.0);
244
tau_ui = 1.0/((alpha_ui + beta_ui)*K_Q10);
245
ui_inf = 1.0/(1.0 + exp((V - 99.45)/27.48));
247
alpha_xr = 0.0003*(V + 14.1)/(1.0 - exp((V + 14.1)/-5.0));
248
beta_xr = 7.3898e-05*(V - 3.3328)/(exp((V -3.3328)/5.1237) - 1.0);
249
tau_xr = 1.0/(alpha_xr + beta_xr);
250
xr_inf = 1.0/(1.0 + exp((V + 14.1)/-6.5));
252
alpha_xs = 4e-05*(V - 19.9)/(1.0 - exp((V - 19.9)/-17.0));
253
beta_xs = 3.5e-05*(V - 19.9)/(exp((V - 19.9)/9.0) - 1.0);
254
tau_xs = 0.5/(alpha_xs + beta_xs);
255
xs_inf = pow((1.0 + exp((V - 19.9)/-12.7)),-0.5);
257
tau_d = (1.0 - exp((V + 10.0)/-6.24))/(0.035*(V + 10.0)*(1.0 + exp((V + 10.0)/-6.24)));
258
d_inf = 1.0/(1.0 + exp((V +10.0)/-8.0));
260
tau_f = 9.0/(0.0197*exp(-0.0337*0.0337*(V + 10.0)*(V + 10.0)) + 0.02);
261
f_inf = 1.0/(1.0 + exp((V + 28.0)/6.9));
263
fca_inf = 1.0/(1.0 + Ca_i/0.00035);
265
Fn = 1e-12*Vrel*I_rel - (5e-13/F)*(0.5*I_CaL - 0.2*I_NaCa);
266
u_inf = 1.0/(1.0 + exp((Fn - 3.4175e-13)/-13.67e-16));
268
tau_v = 1.91 + 2.09/(1.0 + exp((Fn - 3.4175e-13)/-13.67e-16));
269
v_inf = 1.0 - 1.0/(1.0 + exp((Fn - 6.835e-14)/-13.67e-16));
271
tau_w = 6.0*(1.0 - exp((V - 7.9)/-5.0))/((1.0 + 0.3*exp((V - 7.9)/-5.0))*(V - 7.9));
272
w_inf = 1.0 - 1.0/(1.0 + exp((V - 40.0)/-17.0));
275
bool update(const uBlasVector& u, real t, bool end)
285
real m, h, j, oa, oi, ua, ui, xr, xs, d, ff, fca, uu, v, w, Na_i, Ca_i;
286
real Ca_rel, Ca_up, K_i, V;
288
// Ionic currents and gating variables
289
real alpha_m, beta_m, tau_m, m_inf, alpha_h, beta_h, tau_h, h_inf;
290
real alpha_j, beta_j, tau_j, j_inf, alpha_oa, beta_oa, tau_oa, oa_inf;
291
real alpha_oi, beta_oi, tau_oi, oi_inf, alpha_ua, beta_ua, tau_ua, ua_inf;
292
real alpha_ui, beta_ui, tau_ui, ui_inf, alpha_xr, beta_xr, tau_xr, xr_inf;
293
real alpha_xs, beta_xs, tau_xs, xs_inf, tau_d, d_inf, tau_f, f_inf;
294
real fca_inf, u_inf, tau_v, v_inf, tau_w, w_inf, B1, B2;
297
real I_rel, I_CaL, I_NaCa, Fn, sigma, f_NaK, I_NaK, E_Na, I_bNa;
298
real I_Na, I_pCa, E_Ca, I_bCa, I_upleak, I_up, I_tr, E_K, I_K1;
299
real I_to, g_Kur, I_Kur, I_Kr, I_Ks, I_ion;
302
real Cm, R, T, F, z_Na, z_K, z_Ca, Na_o, K_o, Ca_o, K_Q10, tau_fca;
303
real k_rel, g_CaL, gamma, I_NaCamax, K_mNa, K_mCa, k_sat, Vrel, tau_u;
304
real Vi, I_NaKmax, K_mNai, K_mKo, g_bNa, g_Na, Vup, I_pCamax, g_bCa;
305
real I_upmax, K_up, Ca_upmax, Trpn_max, K_mTrpn, Cmdn_max, K_mCmdn;
306
real tau_tr, Csqn_max, K_mCsqn, g_K1, g_to, g_Kr, g_Ks;
307
real Na_e, Ca_e, K_e;
312
// Number of function evaluations
313
unsigned int num_fevals;
322
dolfin_set("ODE tolerance", 10.0);
323
dolfin_set("ODE maximum time step", 100.0);
324
dolfin_set("ODE nonlinear solver", "newton");
325
dolfin_set("ODE linear solver", "iterative");
326
dolfin_set("ODE initial time step", 0.25);
328
//dolfin_set("ODE save solution", false);