LifeV
IonicTenTusscher06.cpp
Go to the documentation of this file.
1 //@HEADER
2 /*
3 *******************************************************************************
4 
5  Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
6  Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University
7 
8  This file is part of LifeV.
9 
10  LifeV is free software; you can redistribute it and/or modify
11  it under the terms of the GNU Lesser General Public License as published by
12  the Free Software Foundation, either version 3 of the License, or
13  (at your option) any later version.
14 
15  LifeV is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public License
21  along with LifeV. If not, see <http://www.gnu.org/licenses/>.
22 
23 *******************************************************************************
24 */
25 //@HEADER
26 
27 /*!
28  @file
29  @brief Ionic model of ten Tusscher-Panfilov 2006
30  @date 01-2013
31  @author Simone Rossi <simone.rossi@epfl.ch>
32 
33  @contributors
34  @mantainer Simone Rossi <simone.rossi@epfl.ch>
35  @last update 01-2013
36  */
37 
38 
39 #include <lifev/electrophysiology/solver/IonicModels/IonicTenTusscher06.hpp>
40 
41 
42 namespace LifeV
43 {
44 
45 // ===================================================
46 //! Constructors
47 // ===================================================
48 IonicTenTusscher06::IonicTenTusscher06() :
49  super ( 19 , 18 ),
50  flag ( MCell )
51 {
52  //Membrane capacitance
53  M_membraneCapacitance = 2.0;// In the paper they give it as 2 but in their code I cannot find it
54 
55  //External concentrations
56  Ko = 5.4;
57  Cao = 2.0;
58  Nao = 140.0;
59 
60  //Intracellular volumes
61  Vc = 0.016404;
62  Vsr = 0.001094;
63  Vss = 0.00005468;
64 
65  //Calcium buffering dynamics
66  Bufc = 0.2;
67  Kbufc = 0.001;
68  Bufsr = 10.;
69  Kbufsr = 0.3;
70  Bufss = 0.4;
71  Kbufss = 0.00025;
72 
73  //Intracellular calcium flux dynamics
74  Vmaxup = 0.006375;
75  Kup = 0.00025;
76  Vrel = 0.102; //40.8;
77  k1_ = 0.15;
78  k2_ = 0.045;
79  k3 = 0.060;
80  k4 = 0.005; //0.000015;
81  EC = 1.5;
82  maxsr = 2.5;
83  minsr = 1.;
84  Vleak = 0.00036;
85  Vxfer = 0.0038;
86 
87 
88 
89  //Constants
90  R = 8314.472;
91  F = 96485.3415;
92  T = 310.0;
93  RTONF = (R * T) / F;
94 
95  //Cellular capacitance
96  M_cellularCapacitance = 0.185;
97 
98  //Parameters for currents
99  //Parameters for IKr
100  Gkr = 0.153;
101  //Parameters for Iks
102  pKNa = 0.03;
103 
104  if (flag == MCell)
105  {
106  Gks = 0.098;
107  }
108  else
109  {
110  Gks = 0.932;
111  }
112  //Parameters for Ik1
113  GK1 = 5.405;
114  //Parameters for Ito
115 
116  if (flag == Endo)
117  {
118  Gks = 0.073;
119  }
120  else
121  {
122  Gto = 0.294;
123  }
124 
125  //Parameters for INa
126  GNa = 14.838;
127  //Parameters for IbNa
128  GbNa = 0.00029;
129  //Parameters for INaK
130  KmK = 1.0;
131  KmNa = 40.0;
132  knak = 2.724;
133  //Parameters for ICaL
134  GCaL = 0.00003980;
135  //Parameters for IbCa
136  GbCa = 0.000592;
137  //Parameters for INaCa
138  knaca = 1000;
139  KmNai = 87.5;
140  KmCa = 1.38;
141  ksat = 0.1;
142  n = 0.35;
143  //Parameters for IpCa
144  GpCa = 0.1238;
145  KpCa = 0.0005;
146  //Parameters for IpK;
147  GpK = 0.0146;
148 
149  inverseVcF2 = 1. / (2.*Vc * F);
150  inverseVcF = 1. / (Vc * F);
151  inversevssF2 = 1. / (2.*Vss * F);
152 
153  //V
154  M_restingConditions.at (0) = -86.2;
155  //m
156  M_restingConditions.at (1) = 0.0;
157  //h
158  M_restingConditions.at (2) = 0.75;
159  //j
160  M_restingConditions.at (3) = 0.75;
161  //d
162  M_restingConditions.at (4) = 0.;
163  //f
164  M_restingConditions.at (5) = 1.;
165  //f2
166  M_restingConditions.at (6) = 1;
167  //fCass
168  M_restingConditions.at (7) = 1.0;
169  //r
170  M_restingConditions.at (8) = 0.;
171  //s
172  M_restingConditions.at (9) = 1.0;
173  //Xr1
174  M_restingConditions.at (10) = 0.0;
175  //Xr2
176  M_restingConditions.at (11) = 1.0;
177  //Xs
178  M_restingConditions.at (12) = 0.0;
179  //Nai
180  M_restingConditions.at (13) = 7.67;
181  //Ki
182  M_restingConditions.at (14) = 138.3;
183  //Cai
184  M_restingConditions.at (15) = 0.00007;
185  //Cass
186  M_restingConditions.at (16) = 0.00007;
187  //Casr
188  M_restingConditions.at (17) = 1.3;
189  //Rprime
190  M_restingConditions.at (18) = 1.0;
191 
192  // M_restingConditions.at(0) = -86.2;
193  // //m
194  // M_restingConditions.at(1) = 0.00165;
195  // //h
196  // M_restingConditions.at(2) = 0.749;
197  // //j
198  // M_restingConditions.at(3) = 0.6788;
199  // //d
200  // M_restingConditions.at(4) = 3.288e-5;
201  // //f
202  // M_restingConditions.at(5) = 0.7026;
203  // //f2
204  // M_restingConditions.at(6) = 0.9526;
205  // //fCass
206  // M_restingConditions.at(7) = 1.0;
207  // //r
208  // M_restingConditions.at(8) = 2.347e-8;
209  // //s
210  // M_restingConditions.at(9) = 0.999998;
211  // //Xr1
212  // M_restingConditions.at(10) = 0.0165;
213  // //Xr2
214  // M_restingConditions.at(11) = 0.473;
215  // //Xs
216  // M_restingConditions.at(12) = 0.0174;
217  // //Nai
218  // M_restingConditions.at(13) = 7.67;
219  // //Ki
220  // M_restingConditions.at(14) = 138.3;
221  // //Cai
222  // M_restingConditions.at(15) = 0.00007;
223  // //Cass
224  // M_restingConditions.at(16) = 0.00007;
225  // //Casr
226  // M_restingConditions.at(17) = 1.3;
227  // //Rprime
228  // M_restingConditions.at(18) = 0.8978;
229 
230 
231 }
232 
233 IonicTenTusscher06::IonicTenTusscher06 ( Teuchos::ParameterList& parameterList ) :
234  super ( 19, 12 )
235 {
236  knak = parameterList.get ("knak", 0.0 );
237  KmNa = parameterList.get ("KmNa", 0.0 );
238  KmK = parameterList.get ("KmK", 0.0 );
239  knaca = parameterList.get ("knaca", 0.0 );
240  KmNai = parameterList.get ("KmNai", 0.0 );
241  KmCa = parameterList.get ("KmCa", 0.0 );
242  ksat = parameterList.get ("ksat", 0.0 );
243  n = parameterList.get ("n", 0.0 );
244 
245 
246  Ko = parameterList.get ("Ko", 0.0 );
247  Cao = parameterList.get ("Cao", 0.0 );
248  Nao = parameterList.get ("Nao", 0.0 );
249 
250 
251  Bufc = parameterList.get ("Bufc", 0.0 );
252  Kbufc = parameterList.get ("Kbufc", 0.0 );
253  Bufsr = parameterList.get ("Bufsr", 0.0 );
254  Kbufsr = parameterList.get ("Kbufsr", 0.0 );
255  Bufss = parameterList.get ("Bufss", 0.0 );
256  Kbufss = parameterList.get ("Kbufss", 0.0 );
257 
258  Vmaxup = parameterList.get ("Vmaxup", 0.0 );
259  Kup = parameterList.get ("Kup", 0.0 );
260  Vrel = parameterList.get ("Vrel", 0.0 );
261  k1_ = parameterList.get ("k1", 0.0 );
262  k2_ = parameterList.get ("k2", 0.0 );
263  k3 = parameterList.get ("k3", 0.0 );
264  k4 = parameterList.get ("k4", 0.0 );
265  EC = parameterList.get ("EC", 0.0 );
266  maxsr = parameterList.get ("maxsr", 0.0 );
267  minsr = parameterList.get ("minsr", 0.0 );
268  Vleak = parameterList.get ("Vleak", 0.0 );
269  Vxfer = parameterList.get ("Vxfer", 0.0 );
270 
271 
272  pKNa = parameterList.get ("pKNa", 0.0 );
273 
274 
275  M_cellularCapacitance = parameterList.get ("Cm", 0.0 );
276  F = parameterList.get ("F", 0.0 );
277  R = parameterList.get ("R", 0.0 );
278  T = parameterList.get ("T", 0.0 );
279 
280 
281  Gkr = parameterList.get ("Gkr", 0.0 );
282  Gks = parameterList.get ("Gks", 0.0 );
283  GK1 = parameterList.get ("GK1", 0.0 );
284  Gto = parameterList.get ("Gto", 0.0 );
285  GNa = parameterList.get ("GNa", 0.0 );
286  GbNa = parameterList.get ("GbNa", 0.0 );
287  GCaL = parameterList.get ("GCal", 0.0 );
288  GbCa = parameterList.get ("GbCa", 0.0 );
289  GpCa = parameterList.get ("GpCa", 0.0 );
290  KpCa = parameterList.get ("KpCa", 0.0 );
291  GpK = parameterList.get ("GpK", 0.0 );
292 
293 
294  Vc = parameterList.get ("Vc", 0.0 );
295  Vsr = parameterList.get ("Vsr", 0.0 );
296  Vss = parameterList.get ("Vss", 0.0 );
297 
298  computeRTONF();
299  computeinverseVcF2();
300  computeinverseVcF();
301  computeinversevssF2();
302 
303  std::map< std::string, WallFlag > WallFlagMap;
304  WallFlagMap["Endo"] = Endo;
305  WallFlagMap["Epi"] = Epi;
306  WallFlagMap["MCell"] = MCell;
307 
308  flag = WallFlagMap[ parameterList.get ( "WallFlag", "MCell" ) ];
309  //V
310  M_restingConditions.at (0) = parameterList.get ("V", -85.423);
311  //m
312  M_restingConditions.at (1) = parameterList.get ("m", 0.00165);
313  //h
314  M_restingConditions.at (2) = parameterList.get ("h", 0.749);
315  //j
316  M_restingConditions.at (3) = parameterList.get ("j", 0.6788);
317  //d
318  M_restingConditions.at (4) = parameterList.get ("d", 3.288e-5);
319  //f
320  M_restingConditions.at (5) = parameterList.get ("f", 0.7026);
321  //f2
322  M_restingConditions.at (6) = parameterList.get ("f2", 0.9526);
323  //fCass
324  M_restingConditions.at (7) = parameterList.get ("fCass", 0.9942);
325  //r
326  M_restingConditions.at (8) = parameterList.get ("r", 2.347e-8);
327  //s
328  M_restingConditions.at (9) = parameterList.get ("s", 0.999998);
329  //Xr1
330  M_restingConditions.at (10) = parameterList.get ("Xr1", 0.0165);
331  //Xr2
332  M_restingConditions.at (11) = parameterList.get ("Xr2", 0.473);
333  //Xs
334  M_restingConditions.at (12) = parameterList.get ("Xs", 0.0174);
335  //Nai
336  M_restingConditions.at (13) = parameterList.get ("Nai", 10.132);
337  //Ki
338  M_restingConditions.at (14) = parameterList.get ("Ki", 138.52);
339  //Cai
340  M_restingConditions.at (15) = parameterList.get ("Cai", 0.000153);
341  //Cass
342  M_restingConditions.at (16) = parameterList.get ("Cass", -85.423);
343  //Casr
344  M_restingConditions.at (17) = parameterList.get ("Casr", 4.272);
345  //Rprime
346  M_restingConditions.at (18) = parameterList.get ("Rprime", 0.8978);
347 }
348 
349 IonicTenTusscher06::IonicTenTusscher06 ( const IonicTenTusscher06& model )
350 {
351  knak = model.knak;
352  KmNa = model.KmNa;
353  KmK = model.KmK;
354  knaca = model.knaca;
355  KmNai = model.KmNai;
356  KmCa = model.KmCa;
357  ksat = model.ksat;
358  n = model.n;
359 
360 
361  Ko = model.Ko;
362  Cao = model.Cao;
363  Nao = model.Nao;
364 
365 
366  Bufc = model.Bufc;
367  Kbufc = model.Kbufc;
368  Bufsr = model.Bufsr;
369  Kbufsr = model.Kbufsr;
370  Bufss = model.Bufss;
371  Kbufss = model.Kbufss;
372 
373  Vmaxup = model.Vmaxup;
374  Kup = model.Kup;
375  Vrel = model.Vrel;
376  k1_ = model.k1_;
377  k2_ = model.k2_;
378  k3 = model.k3;
379  k4 = model.k4;
380  EC = model.EC;
381  maxsr = model.maxsr;
382  minsr = model.minsr;
383  Vleak = model.Vleak;
384  Vxfer = model.Vxfer;
385 
386 
387  pKNa = model.pKNa;
388 
389  F = model.F;
390  R = model.R;
391  T = model.T;
392 
393 
394  Gkr = model.Gkr;
395  Gks = model.Gks;
396  GK1 = model.GK1;
397  Gto = model.Gto;
398  GNa = model.GNa;
399  GbNa = model.GbNa;
400  GCaL = model.GCaL;
401  GbCa = model.GbCa;
402  GpCa = model.GpCa;
403  KpCa = model.KpCa;
404  GpK = model.GpK;
405 
406 
407  Vc = model.Vc;
408  Vsr = model.Vsr;
409  Vss = model.Vss;
410 
411  computeRTONF();
412  computeinverseVcF2();
413  computeinverseVcF();
414  computeinversevssF2();
415 
416  flag = model.flag;
417 
418  M_numberOfEquations = model.M_numberOfEquations;
419  M_numberOfGatingVariables = model.M_numberOfGatingVariables;
420  M_restingConditions = model.M_restingConditions;
421 
422 }
423 
424 // ===================================================
425 //! Operator
426 // ===================================================
427 IonicTenTusscher06& IonicTenTusscher06::operator= ( const IonicTenTusscher06& model )
428 {
429  knak = model.knak;
430  KmNa = model.KmNa;
431  KmK = model.KmK;
432  knaca = model.knaca;
433  KmNai = model.KmNai;
434  KmCa = model.KmCa;
435  ksat = model.ksat;
436  n = model.n;
437 
438 
439  Ko = model.Ko;
440  Cao = model.Cao;
441  Nao = model.Nao;
442 
443 
444  Bufc = model.Bufc;
445  Kbufc = model.Kbufc;
446  Bufsr = model.Bufsr;
447  Kbufsr = model.Kbufsr;
448  Bufss = model.Bufss;
449  Kbufss = model.Kbufss;
450 
451  Vmaxup = model.Vmaxup;
452  Kup = model.Kup;
453  Vrel = model.Vrel;
454  k1_ = model.k1_;
455  k2_ = model.k2_;
456  k3 = model.k3;
457  k4 = model.k4;
458  EC = model.EC;
459  maxsr = model.maxsr;
460  minsr = model.minsr;
461  Vleak = model.Vleak;
462  Vxfer = model.Vxfer;
463 
464 
465  pKNa = model.pKNa;
466 
467  F = model.F;
468  R = model.R;
469  T = model.T;
470 
471 
472  Gkr = model.Gkr;
473  Gks = model.Gks;
474  GK1 = model.GK1;
475  Gto = model.Gto;
476  GNa = model.GNa;
477  GbNa = model.GbNa;
478  GCaL = model.GCaL;
479  GbCa = model.GbCa;
480  GpCa = model.GpCa;
481  KpCa = model.KpCa;
482  GpK = model.GpK;
483 
484 
485  Vc = model.Vc;
486  Vsr = model.Vsr;
487  Vss = model.Vss;
488 
489  computeRTONF();
490  computeinverseVcF2();
491  computeinverseVcF();
492  computeinversevssF2();
493 
494  flag = model.flag;
495 
496  M_numberOfEquations = model.M_numberOfEquations;
497  M_numberOfGatingVariables = model.M_numberOfGatingVariables;
498  M_restingConditions = model.M_restingConditions;
499 
500  return *this;
501 }
502 
503 
504 // ===================================================
505 //! Methods
506 // ===================================================
507 void IonicTenTusscher06::computeGatingRhs ( const std::vector<Real>& v,
508  std::vector<Real>& rhs )
509 {
510  Real V = v[0];
511  Real m = v[1];
512  Real h = v[2];
513  Real j = v[3];
514  Real d = v[4];
515  Real f = v[5];
516  Real f2 = v[6];
517  Real fcass = v[7];
518  Real r = v[8];
519  Real s = v[9];
520  Real xr1 = v[10];
521  Real xr2 = v[11];
522  Real xs = v[12];
523  Real Nai = v[13];
524  Real Ki = v[14];
525  Real Cai = v[15];
526  Real CaSS = v[16];
527  Real CaSR = v[17];
528  Real RR = v[18];
529 
530  //m
531  rhs[0] = dM (V, m);
532  //h
533  rhs[1] = dH (V, h);
534  //j
535  rhs[2] = dJ (V, j);
536  //d
537  rhs[3] = dD (V, d);
538  //f
539  rhs[4] = dF (V, f);
540  //f2
541  rhs[5] = dF2 (V, f2);
542  //fCass
543  rhs[6] = dFCaSS (V, fcass);
544  //r
545  rhs[7] = dR (V, r);
546  //s
547  rhs[8] = dS (V, s);
548  //Xr1
549  rhs[9] = dXr1 (V, xr1);
550  //Xr2
551  rhs[10] = dXr2 (V, xr2);
552  //Xs
553  rhs[11] = dXs (V, xs);
554  //Nai
555  rhs[12] = dNai (V, m, h, j, Nai, Cai);
556  //Ki
557  rhs[13] = dKi (V, r, s, xr1, xr2, xs, Ki, Nai);
558  //Cai
559  rhs[14] = dCai (V, Nai, Cai, CaSR, CaSS);
560  //Cass
561  rhs[15] = dCaSS (Cai, CaSR, CaSS, RR, V, d, f, f2, fcass);
562  //Casr
563  rhs[16] = dCaSR (Cai, CaSR, CaSS, RR);
564  //Rprime
565  rhs[17] = dRR (CaSR, CaSS, RR);
566 }
567 
568 void IonicTenTusscher06::computeNonGatingRhs ( const std::vector<Real>& /*v*/,
569  std::vector<Real>& /*rhs*/ )
570 {
571 }
572 
573 void IonicTenTusscher06::computeRhs ( const std::vector<Real>& v,
574  std::vector<Real>& rhs )
575 {
576  Real V = v[0];
577  Real m = v[1];
578  Real h = v[2];
579  Real j = v[3];
580  Real d = v[4];
581  Real f = v[5];
582  Real f2 = v[6];
583  Real fcass = v[7];
584  Real r = v[8];
585  Real s = v[9];
586  Real xr1 = v[10];
587  Real xr2 = v[11];
588  Real xs = v[12];
589  Real Nai = v[13];
590  Real Ki = v[14];
591  Real Cai = v[15];
592  Real CaSS = v[16];
593  Real CaSR = v[17];
594  Real RR = v[18];
595 
596  //V
597  rhs[0] = - Itot (V, m, h, j, d, f, f2, fcass, r, s, xr1, xr2, xs, Nai, Ki, Cai, CaSS );
598  //m
599  rhs[1] = dM (V, m);
600  //h
601  rhs[2] = dH (V, h);
602  //j
603  rhs[3] = dJ (V, j);
604  //d
605  rhs[4] = dD (V, d);
606  //f
607  rhs[5] = dF (V, f);
608  //f2
609  rhs[6] = dF2 (V, f2);
610  //fCass
611  rhs[7] = dFCaSS (V, fcass);
612  //r
613  rhs[8] = dR (V, r);
614  //s
615  rhs[9] = dS (V, s);
616  //Xr1
617  rhs[10] = dXr1 (V, xr1);
618  //Xr2
619  rhs[11] = dXr2 (V, xr2);
620  //Xs
621  rhs[12] = dXs (V, xs);
622  //Nai
623  rhs[13] = dNai (V, m, h, j, Nai, Cai);
624  //Ki
625  rhs[14] = dKi (V, r, s, xr1, xr2, xs, Ki, Nai);
626  //Cai
627  rhs[15] = dCai (V, Nai, Cai, CaSR, CaSS);
628  //Cass
629  rhs[16] = dCaSS (Cai, CaSR, CaSS, RR, V, d, f, f2, fcass);
630  //Casr
631  rhs[17] = dCaSR (Cai, CaSR, CaSS, RR);
632  //Rprime
633  rhs[18] = dRR (CaSR, CaSS, RR);
634 
635 }
636 
637 
638 Real IonicTenTusscher06::computeLocalPotentialRhs ( const std::vector<Real>& v)
639 {
640  Real dPotential (0.0);
641 
642  Real V = v[0];
643  Real m = v[1];
644  Real h = v[2];
645  Real j = v[3];
646  Real d = v[4];
647  Real f = v[5];
648  Real f2 = v[6];
649  Real fCass = v[7];
650  Real r = v[8];
651  Real s = v[9];
652  Real Xr1 = v[10];
653  Real Xr2 = v[11];
654  Real Xs = v[12];
655  Real Nai = v[13];
656  Real Ki = v[14];
657  Real Cai = v[15];
658  Real Cass = v[16];
659 
660  dPotential = -Itot (V, m, h, j, d, f, f2, fCass, r, s, Xr1, Xr2, Xs, Nai, Ki, Cai, Cass );
661 
662  return dPotential;
663 }
664 
665 
666 void IonicTenTusscher06::computeGatingVariablesWithRushLarsen ( std::vector<Real>& v, const Real dt )
667 {
668  Real V = v[0];
669  Real m = v[1];
670  Real h = v[2];
671  Real j = v[3];
672  Real d = v[4];
673  Real f = v[5];
674  Real f2 = v[6];
675  Real fcass = v[7];
676  Real r = v[8];
677  Real s = v[9];
678  Real xr1 = v[10];
679  Real xr2 = v[11];
680  Real xs = v[12];
681  Real Nai = v[13];
682  Real Ki = v[14];
683  Real Cai = v[15];
684  Real CaSS = v[16];
685  Real CaSR = v[17];
686  Real RR = v[18];
687 
688  v[1] = M_INF (V) - ( M_INF (V) - m ) * std::exp (- dt / TAU_M (V) );
689  v[2] = H_INF (V) - ( H_INF (V) - h ) * std::exp (- dt / TAU_H (V) );
690  v[3] = J_INF (V) - ( J_INF (V) - j ) * std::exp (- dt / TAU_J (V) );
691  v[4] = D_INF (V) - ( D_INF (V) - d ) * std::exp (- dt / TAU_D (V) );
692  v[5] = F_INF (V) - ( F_INF (V) - f ) * std::exp (- dt / TAU_F (V) );
693  v[6] = F2_INF (V) - ( F2_INF (V) - f2 ) * std::exp (- dt / TAU_F2 (V) );
694  v[7] = FCaSS_INF (CaSS) - ( FCaSS_INF (CaSS) - fcass ) * std::exp ( -dt / TAU_FCaSS (CaSS) );
695  v[8] = R_INF (V) - ( R_INF (V) - r ) * std::exp (- dt / TAU_R (V) );
696  v[9] = S_INF (V) - ( S_INF (V) - s ) * std::exp (- dt / TAU_S (V) );
697  v[10] = Xr1_INF (V) - ( Xr1_INF (V) - xr1 ) * std::exp (- dt / TAU_Xr1 (V) );
698  v[11] = Xr2_INF (V) - ( Xr2_INF (V) - xr2 ) * std::exp (- dt / TAU_Xr2 (V) );
699  v[12] = Xs_INF (V) - ( Xs_INF (V) - xs ) * std::exp (- dt / TAU_Xs (V) );
700  v[13] = solveNai (V, m, h, j, Nai, Cai, dt);
701  v[14] = solveKi (V, r, s, xr1, xr2, xs, Nai, Ki, dt);
702  v[15] = solveCai (V, Nai, Cai, CaSR, CaSS, dt);
703  v[16] = solveCaSS (Cai, CaSR, CaSS, RR, V, d, f, f2, fcass, dt);
704  v[17] = solveCaSR (Cai, CaSR, CaSS, RR, dt);
705  v[18] = solveRR (CaSR, CaSS, RR, dt);
706 
707 
708 }
709 
710 void IonicTenTusscher06::showMe()
711 {
712  std::cout << "\n\n************************************";
713  std::cout << "\n\tHi, I'm the Ten Tusscher model";
714  std::cout << "\n\t I've so many parameters that I don't think it's a good idea to display them all\n\n";
715  std::cout << "\n\tPlease use a getter, or implement this method otherwise.";
716  std::cout << "\n************************************\n\n";
717 }
718 
719 void IonicTenTusscher06::solveOneStep (std::vector<Real>& v, Real dt)
720 {
721  Real V = v[0];
722  Real m = v[1];
723  Real h = v[2];
724  Real j = v[3];
725  Real d = v[4];
726  Real f = v[5];
727  Real f2 = v[6];
728  Real fcass = v[7];
729  Real r = v[8];
730  Real s = v[9];
731  Real xr1 = v[10];
732  Real xr2 = v[11];
733  Real xs = v[12];
734  Real Nai = v[13];
735  Real Ki = v[14];
736  Real Cai = v[15];
737  Real CaSS = v[16];
738  Real CaSR = v[17];
739  Real RR = v[18];
740 
741  computeGatingVariablesWithRushLarsen (v, dt);
742 
743  v[0] = solveV (V, m, h, j, d, f, f2, fcass, r, s, xr1, xr2, xs, Nai, Ki, Cai, CaSS, dt);
744  v[13] = solveNai (V, m, h, j, Nai, Cai, dt);
745  v[14] = solveKi (V, r, s, xr1, xr2, xs, Nai, Ki, dt);
746  v[15] = solveCai (V, Nai, Cai, CaSR, CaSS, dt);
747  v[16] = solveCaSS (Cai, CaSR, CaSS, RR, V, d, f, f2, fcass, dt);
748  v[17] = solveCaSR (Cai, CaSR, CaSS, RR, dt);
749  v[18] = solveRR (CaSR, CaSS, RR, dt);
750 
751 
752 
753 
754 }
755 
756 
757 void IonicTenTusscher06::showCurrents (std::vector<Real>& v)
758 {
759  Real V = v[0];
760  Real m = v[1];
761  Real h = v[2];
762  Real j = v[3];
763  Real d = v[4];
764  Real f = v[5];
765  Real f2 = v[6];
766  Real fcass = v[7];
767  Real r = v[8];
768  Real s = v[9];
769  Real xr1 = v[10];
770  Real xr2 = v[11];
771  Real xs = v[12];
772  Real Nai = v[13];
773  Real Ki = v[14];
774  Real Cai = v[15];
775  Real CaSS = v[16];
776 
777  showCurrents ( V, m, h, j, d, f, f2, fcass, r, s, xr1, xr2, xs, Nai, Ki, Cai, CaSS);
778 }
779 
780 void IonicTenTusscher06::showCurrents (Real V, Real m, Real h, Real j, Real d, Real f, Real f2, Real fcass,
781  Real r, Real s, Real xr1, Real xr2, Real xs, Real Nai, Real Ki,
782  Real Cai, Real CaSS)
783 {
784  std::cout << "\nIKr = " << IKr (V, xr1, xr2, Ki);
785  std::cout << "\nIKs = " << IKs (V, xs, Ki, Nai);
786  std::cout << "\nIK1 = " << IK1 (V, Ki);
787  std::cout << "\nIto = " << Ito (V, r, s, Ki);
788  std::cout << "\nINa = " << INa (V, m, h, j, Nai);
789  std::cout << "\nIbNa = " << IbNa (V, Nai);
790  std::cout << "\nICaL = " << ICaL (V, d, f, f2, fcass, CaSS);
791  std::cout << "\nIbCa = " << IbCa (V, Cai);
792  std::cout << "\nINaK = " << INaK (V, Nai);
793  std::cout << "\nINaCa = " << INaCa (V, Nai, Cai);
794  std::cout << "\nIpCa = " << IpCa (Cai);
795  std::cout << "\nIKpK = " << IpK (V, Ki);
796  std::cout << "\n";
797 }
798 
799 
800 }
ETCurrentFE::updateInverseJacobian
void updateInverseJacobian(const UInt &iQuadPt)
Definition: ETCurrentFE.cpp:405
LifeV::Real
double Real
Generic real data.
Definition: LifeV.hpp:175