← Back to branch summary

~njansson/dolfin/hpc

« back to all changes in this revision

Viewing changes to src/demo/ode/courtemanche/main.cpp

  • Committer: Johannes Ring
  • Date: 2008-03-05 22:43:06 UTC
  • Revision ID: johannr@simula.no-20080305224306-2npsdyhfdpl2esji
The BIG commit!

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/demo/ode/courtemanche/main.cpp
1
 
// Copyright (C) Glenn Terje Lines, Ola Skavhaug and Simula Research Laboratory.
2
 
// Licensed under the GNU LGPL Version 2.1.
3
 
//
4
 
// Original code copied from PyCC.
5
 
// Modified by Anders Logg 2006.
6
 
//
7
 
// First added:  2006-05-24
8
 
// Last changed: 2006-08-21
9
 
//
10
 
// This demo solves the Courtemanche model for cardiac excitation.
11
 
 
12
 
#include <dolfin.h>
13
 
 
14
 
using namespace dolfin;
15
 
 
16
 
class Courtemanche : public ODE
17
 
{
18
 
public:
19
 
  
20
 
  Courtemanche() : ODE(21, 300.0)
21
 
  {
22
 
    // Set parameters
23
 
    Cm        = 100.0;
24
 
    R         = 8.3143; 
25
 
    T         = 310; 
26
 
    F         = 96.4867; 
27
 
    z_Na      = 1.0; 
28
 
    z_K       = 1.0; 
29
 
    z_Ca      = 2.0;
30
 
    Na_o      = 140.0; 
31
 
    K_o       = 5.4; 
32
 
    Ca_o      = 1.8; 
33
 
    K_Q10     = 3.0;
34
 
    tau_fca   = 2.0;
35
 
    k_rel     = 30.0;
36
 
    g_CaL     = 0.12375;
37
 
    gamma     = 0.35;
38
 
    I_NaCamax = 1600.0;
39
 
    K_mNa     = 87.5;
40
 
    K_mCa     = 1.38;
41
 
    k_sat     = 0.1;
42
 
    Vrel      = 96.48;
43
 
    tau_u     = 8.0;
44
 
    Vi        = 13668.0;
45
 
    I_NaKmax  = 0.59933874;
46
 
    K_mNai    = 10.0;
47
 
    //K_o=5.4;
48
 
    K_mKo     = 1.5;
49
 
    g_bNa     = 0.0006744375;
50
 
    g_Na      = 7.8;
51
 
    Vup       = 1109.52;
52
 
    I_pCamax  = 0.275;
53
 
    g_bCa     = 0.001131;
54
 
    I_upmax   = 0.005;
55
 
    K_up      = 0.00092;
56
 
    Ca_upmax  = 15.0;
57
 
    Trpn_max  = 0.070;
58
 
    K_mTrpn   = 0.0005;
59
 
    Cmdn_max  = 0.050;
60
 
    K_mCmdn   = 0.00238;
61
 
    tau_tr    = 180.0;
62
 
    Csqn_max  = 10.0;
63
 
    K_mCsqn   = 0.8;
64
 
    g_K1      = 0.09;
65
 
    g_to      = 0.1652;
66
 
    g_Kr      = 0.029411765;
67
 
    g_Ks      = 0.12941176;
68
 
    g_Na      = 7.8;
69
 
    ist       = 0.0;
70
 
 
71
 
    num_fevals = 0;
72
 
    VT = 0.0;
73
 
  }
74
 
  
75
 
  ~Courtemanche()
76
 
  {
77
 
    message("Function evaluations:  %d", num_fevals);
78
 
    message("Potential at end time: %.6f", VT);
79
 
  }
80
 
 
81
 
  void u0(uBlasVector& u)
82
 
  {
83
 
    // Set initial data
84
 
    u(0)  = -85.0; 
85
 
    u(1)  = 2.91e-3; 
86
 
    u(2)  = 9.65e-1;
87
 
    u(3)  = 9.78e-1;
88
 
    u(4)  = 3.04e-2;
89
 
    u(5)  = 9.99e-1;
90
 
    u(6)  = 4.96e-3;
91
 
    u(7)  = 9.99e-1;
92
 
    u(8)  = 3.29e-5;
93
 
    u(9)  = 1.87e-2;
94
 
    u(10) = 1.37e-4;
95
 
    u(11) = 9.99e-1; 
96
 
    u(12) = 7.75e-1;
97
 
    u(13) = 0.0;
98
 
    u(14) = 1.0;
99
 
    u(15) = 9.99e-1;
100
 
    u(16) = 11.2;
101
 
    u(17) = 1.02e-4;
102
 
    u(18) = 1.49;
103
 
    u(19) = 1.49;
104
 
    u(20) = 139.0;
105
 
 
106
 
    // Initial kick
107
 
    u(0) = -25.0;
108
 
  }
109
 
  
110
 
  void f(const uBlasVector& u, real t, uBlasVector& y)
111
 
  {
112
 
    computeCurrents(u);
113
 
    computeGateCoefficients(u);
114
 
 
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);
132
 
    y(17) = B1/B2;
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);
136
 
 
137
 
    num_fevals++;
138
 
  }
139
 
 
140
 
  void computeCurrents(const uBlasVector& u)
141
 
  {
142
 
    V      = u(0);
143
 
    m      = u(1);
144
 
    h      = u(2);
145
 
    j      = u(3);
146
 
    oa     = u(4);
147
 
    oi     = u(5);
148
 
    ua     = u(6);
149
 
    ui     = u(7);
150
 
    xr     = u(8);
151
 
    xs     = u(9);
152
 
    d      = u(10);
153
 
    ff     = u(11);
154
 
    fca    = u(12);
155
 
    uu     = u(13);
156
 
    v      = u(14);
157
 
    w      = u(15);
158
 
    Na_i   = u(16);
159
 
    Ca_i   = u(17);
160
 
    Ca_rel = u(18);
161
 
    Ca_up  = u(19);
162
 
    K_i    = u(20);
163
 
 
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));
190
 
  }
191
 
  
192
 
  void computeGateCoefficients(const uBlasVector& u)
193
 
  {
194
 
    V = u(0);
195
 
    
196
 
    if ( V == -47.13 )
197
 
        alpha_m = 3.2;
198
 
    else
199
 
        alpha_m = 0.32*(V + 47.13)/(1.0 - exp(-0.1*(V + 47.13)));
200
 
 
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;
204
 
    if (V >= -40.0){
205
 
        alpha_h = 0.0;
206
 
        beta_h  = 1.0/(0.13*(1.0 + exp((V + 10.66)/-11.1)));
207
 
    } else {
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);
210
 
    }
211
 
 
212
 
    tau_h = 1.0/(alpha_h + beta_h);
213
 
    h_inf = alpha_h*tau_h;
214
 
      
215
 
    if ( V >= -40.0 )
216
 
    {
217
 
        alpha_j = 0.0;
218
 
        beta_j  = 0.3*(exp(-2.535e-7*V))/(1.0 + exp(-0.1*(V +32.0)));
219
 
    } else {
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)));
222
 
    }
223
 
 
224
 
    tau_j = 1.0/(alpha_j + beta_j);
225
 
    j_inf = alpha_j*tau_j;
226
 
      
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));
231
 
      
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));
236
 
      
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));
241
 
     
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));
246
 
      
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));
251
 
      
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);
256
 
     
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));
259
 
      
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));
262
 
     
263
 
    fca_inf  = 1.0/(1.0 + Ca_i/0.00035);
264
 
     
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));
267
 
     
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));
270
 
     
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));
273
 
  }
274
 
  
275
 
  bool update(const uBlasVector& u, real t, bool end)
276
 
  {
277
 
    if ( end )
278
 
      VT = u(0);
279
 
    return true;
280
 
  }
281
 
  
282
 
private:
283
 
  
284
 
  // State varibles
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;
287
 
  
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;
295
 
  
296
 
  // Membrane currents
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;
300
 
  
301
 
  // Gate coefficients
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;
308
 
 
309
 
  // Stimulus current
310
 
  real ist;
311
 
  
312
 
  // Number of function evaluations
313
 
  unsigned int num_fevals;
314
 
 
315
 
  // Value at end time
316
 
  real VT;
317
 
 
318
 
};
319
 
 
320
 
int main()
321
 
{
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);
327
 
 
328
 
  //dolfin_set("ODE save solution", false);
329
 
 
330
 
  Courtemanche ode;
331
 
  ode.solve();
332
 
 
333
 
  return 0;
334
 
}