LifeV
IonicFox.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 based on Fox model.
30  @date 04-2013
31  @author Marie Dupraz <dupraz.marie@gmail.com>
32 
33  @contributors
34  @mantainer Marie Dupraz <dupraz.marie@gmail.com>
35  @last update 04-2013
36  */
37 
38 #include <lifev/electrophysiology/solver/IonicModels/IonicFox.hpp>
39 #include <Teuchos_RCP.hpp>
40 #include <Teuchos_ParameterList.hpp>
41 #include "Teuchos_XMLParameterListHelpers.hpp"
42 
43 #include <cmath>
44 
45 
46 namespace LifeV
47 {
48 // ===================================================
49 //! Constructors
50 // ===================================================
51 IonicFox::IonicFox() :
52  super ( 13, 7 ),
53  M_ACap ( 1.534e-4 ),
54  M_VMyo ( 25.84e-6 ),
55  M_Vup ( 0.1 ),
56  M_Vsr ( 2e-6 ),
57  M_NaO ( 138.0 ),
58  M_CaO ( 2000.0 ),
59  M_K0 ( 4.0 ),
60  M_NaIn ( 10.0 ),
61  M_KIn ( 149.4 ),
62  M_CmdnTot ( 10.0 ),
63  M_CsqnTot ( 10000.0 ),
64  M_KmCmdn ( 2.0 ),
65  M_KmCsqn ( 600.0 ),
66  M_Cm ( 1.0 ),
67  M_F ( 96.5 ),
68  M_T ( 310.0 ),
69  M_R ( 8.314 ),
70  M_GNa ( 12.8 ),
71  M_GKp ( 0.002216 ),
72  M_KmNa ( 87.5 ),
73  M_KmCa ( 1380.0 ),
74  M_kSat ( 0.2 ),
75  M_kNaCa ( 1500.0 ),
76  M_eta ( 0.35 ),
77  M_INaK ( 0.693 ),
78  M_KmNai ( 10.0 ),
79  M_KmK0 ( 1.5 ),
80  M_IpCa ( 0.05 ),
81  M_KmPCa ( 0.05 ),
82  M_GCab ( 0.0003842 ),
83  M_GNab ( 0.0031 ),
84  M_KmUp ( 0.32 ),
85  M_PCa ( 0.0000226 ),
86  M_ICaHalf ( -0.265 ),
87  M_GK1 ( 2.8 ),
88  M_GKr ( 0.0136 ),
89  M_GKs ( 0.0245 ),
90  M_Gt0 ( 0.23815 ),
91  M_PCaK ( 5.79e-7 ),
92  M_Prel ( 6.0 ),
93  M_Pleak ( 0.000001 ),
94  M_KmfCa ( 0.18 ),
95  M_KmK1 ( 13.0 )
96 {
97  M_restingConditions.at (0) = - 94.7;
98  M_restingConditions.at (1) = 2.4676e-4;
99  M_restingConditions.at (2) = 0.99869;
100  M_restingConditions.at (3) = 0.99887;
101  M_restingConditions.at (4) = 0.229;
102  M_restingConditions.at (5) = 0.0001;
103  M_restingConditions.at (6) = 3.742e-5;
104  M_restingConditions.at (7) = 1.0;
105  M_restingConditions.at (8) = 0.983;
106  M_restingConditions.at (9) = 0.0001;
107  M_restingConditions.at (10) = 0.942;
108  M_restingConditions.at (11) = 0.0472;
109  M_restingConditions.at (12) = 320.0;
110 }
111 
112 IonicFox::IonicFox ( Teuchos::ParameterList& parameterList ) :
113  super ( 13, 7 )
114 {
115  M_ACap = parameterList.get ( "areaCap", 1.534e-4 );
116  M_VMyo = parameterList.get ( "volMyo", 25.84e-6 );
117  M_Vup = parameterList.get ( "volUp", 0.1 );
118  M_Vsr = parameterList.get ( "volSR", 2e-6 );
119  M_NaO = parameterList.get ( "concNa0", 138.0 );
120  M_CaO = parameterList.get ( "concCa0", 2000.0 );
121  M_K0 = parameterList.get ( "concK0", 4.0 );
122  M_NaIn = parameterList.get ( "concNaIn", 10.0 );
123  M_KIn = parameterList.get ( "concKIn", 149.4 );
124  M_CmdnTot = parameterList.get ( "cmdnTot", 10.0 );
125  M_CsqnTot = parameterList.get ( "csqnTot", 10000.0 );
126  M_KmCmdn = parameterList.get ( "constmCmdn", 2.0 );
127  M_KmCsqn = parameterList.get ( "constmCsqn", 600.0 );
128  M_Cm = parameterList.get ( "capMem", 1.0 );
129  M_F = parameterList.get ( "farad", 96.5 );
130  M_T = parameterList.get ( "temp", 310.0 );
131  M_R = parameterList.get ( "gasConst", 8.314 );
132  M_GNa = parameterList.get ( "maxCondNa", 12.8 );
133  M_GKp = parameterList.get ( "maxCondKp", 0.002216 );
134  M_KmNa = parameterList.get ( "constmNa", 87.5 );
135  M_KmCa = parameterList.get ( "constmCa", 1380.0 );
136  M_kSat = parameterList.get ( "kSat", 0.2 );
137  M_kNaCa = parameterList.get ( "kNaCa", 1500.0 );
138  M_eta = parameterList.get ( "eta", 0.35 );
139  M_INaK = parameterList.get ( "courNaK", 0.693 );
140  M_KmNai = parameterList.get ( "constmNai", 10.0 );
141  M_KmK0 = parameterList.get ( "constmK0", 1.5 );
142  M_IpCa = parameterList.get ( "courpCa", 0.05 );
143  M_KmPCa = parameterList.get ( "constmpCa", 0.05 );
144  M_GCab = parameterList.get ( "maxCondCab", 0.0003842 );
145  M_GNab = parameterList.get ( "maxCondNab", 0.0031 );
146  M_KmUp = parameterList.get ( "constmUp", 0.32 );
147  M_PCa = parameterList.get ( "permCa", 0.0000226 );
148  M_ICaHalf = parameterList.get ( "courCaHalf", -0.265 );
149  M_GK1 = parameterList.get ( "maxCondK1", 2.8 );
150  M_GKr = parameterList.get ( "maxCondKr", 0.0136 );
151  M_GKs = parameterList.get ( "maxCondKs", 0.0245 );
152  M_Gt0 = parameterList.get ( "maxCondt0", 0.23815 );
153  M_PCaK = parameterList.get ( "permCaK", 5.79e-7 );
154  M_Prel = parameterList.get ( "permrel", 6.0 );
155  M_Pleak = parameterList.get ( "permleak", 0.000001 );
156  M_KmfCa = parameterList.get ( "constmfCa", 0.18 );
157  M_KmK1 = parameterList.get ( "constmK1", 13.0 );
158 }
159 
160 IonicFox::IonicFox ( const IonicFox& model )
161 {
162  M_ACap = model.M_ACap;
163  M_VMyo = model.M_VMyo;
164  M_Vup = model.M_Vup;
165  M_Vsr = model.M_Vsr;
166  M_NaO = model.M_NaO;
167  M_CaO = model.M_CaO;
168  M_K0 = model.M_K0;
169  M_NaIn = model.M_NaIn;
170  M_KIn = model.M_KIn;
171  M_CmdnTot = model.M_CmdnTot;
172  M_CsqnTot = model.M_CsqnTot;
173  M_KmCmdn = model.M_KmCmdn;
174  M_KmCsqn = model.M_KmCsqn;
175  M_Cm = model.M_Cm;
176  M_F = model.M_F;
177  M_T = model.M_T;
178  M_R = model.M_R;
179  M_GNa = model.M_GNa;
180  M_GKp = model.M_GKp;
181  M_KmNa = model.M_KmNa;
182  M_KmCa = model.M_KmCa;
183  M_kSat = model.M_kSat;
184  M_kNaCa = model.M_kNaCa;
185  M_eta = model.M_eta;
186  M_INaK = model.M_INaK;
187  M_KmNai = model.M_KmNai;
188  M_KmK0 = model.M_KmK0;
189  M_IpCa = model.M_IpCa;
190  M_KmPCa = model.M_KmPCa;
191  M_GCab = model.M_GCab;
192  M_GNab = model.M_GNab;
193  M_KmUp = model.M_KmUp;
194  M_PCa = model.M_PCa;
195  M_ICaHalf = model.M_ICaHalf;
196  M_GK1 = model.M_GK1;
197  M_GKr = model.M_GKr;
198  M_GKs = model.M_GKs;
199  M_Gt0 = model.M_Gt0;
200  M_PCaK = model.M_PCaK;
201  M_Prel = model.M_Prel;
202  M_Pleak = model.M_Pleak;
203  M_KmfCa = model. M_KmfCa;
204  M_KmK1 = model. M_KmK1;
205 
206  M_numberOfEquations = model.M_numberOfEquations;
207  M_restingConditions = model.M_restingConditions;
208  M_numberOfGatingVariables = model.M_numberOfGatingVariables;
209 }
210 
211 // ===================================================
212 //! Operator
213 // ===================================================
214 IonicFox& IonicFox::operator= ( const IonicFox& model )
215 {
216  M_ACap = model.M_ACap;
217  M_VMyo = model.M_VMyo;
218  M_Vup = model.M_Vup;
219  M_Vsr = model.M_Vsr;
220  M_NaO = model.M_NaO;
221  M_CaO = model.M_CaO;
222  M_K0 = model.M_K0;
223  M_NaIn = model.M_NaIn;
224  M_KIn = model.M_KIn;
225  M_CmdnTot = model.M_CmdnTot;
226  M_CsqnTot = model.M_CsqnTot;
227  M_KmCmdn = model.M_KmCmdn;
228  M_KmCsqn = model.M_KmCsqn;
229  M_Cm = model.M_Cm;
230  M_F = model.M_F;
231  M_T = model.M_T;
232  M_R = model.M_R;
233  M_GNa = model.M_GNa;
234  M_GKp = model.M_GKp;
235  M_KmNa = model.M_KmNa;
236  M_KmCa = model.M_KmCa;
237  M_kSat = model.M_kSat;
238  M_kNaCa = model.M_kNaCa;
239  M_eta = model.M_eta;
240  M_INaK = model.M_INaK;
241  M_KmNai = model.M_KmNai;
242  M_KmK0 = model.M_KmK0;
243  M_IpCa = model.M_IpCa;
244  M_KmPCa = model.M_KmPCa;
245  M_GCab = model.M_GCab;
246  M_GNab = model.M_GNab;
247  M_KmUp = model.M_KmUp;
248  M_PCa = model.M_PCa;
249  M_ICaHalf = model.M_ICaHalf;
250  M_GK1 = model.M_GK1;
251  M_GKr = model.M_GKr;
252  M_GKs = model.M_GKs;
253  M_Gt0 = model.M_Gt0;
254  M_PCaK = model.M_PCaK;
255  M_Prel = model.M_Prel;
256  M_Pleak = model.M_Pleak;
257  M_KmfCa = model. M_KmfCa;
258  M_KmK1 = model. M_KmK1;
259 
260  M_numberOfEquations = model.M_numberOfEquations;
261  M_restingConditions = model.M_restingConditions;
262  M_numberOfGatingVariables = model.M_numberOfGatingVariables;
263 
264  return *this;
265 }
266 
267 
268 // ===================================================
269 //! Methods
270 // ===================================================
271 //Only gating variables
272 void IonicFox::computeGatingRhs ( const std::vector<Real>& v, std::vector<Real>& rhs )
273 {
274  std::vector<Real> gatingRhs ( computeLocalGatingRhs (v) );
275  std::vector<Real> subSysCaRhs ( computeLocalSubSysCaRhs (v) );
276 
277  std::copy ( gatingRhs.begin(), gatingRhs.end(), rhs.begin() );
278  std::copy ( subSysCaRhs.begin(), subSysCaRhs.end() - 1, rhs.begin() + gatingRhs.size() );
279 }
280 
281 //Potential and gating variables
282 void IonicFox::computeRhs (const std::vector<Real>& v, std::vector<Real>& rhs )
283 {
284  std::vector<Real> gatingRhs ( computeLocalGatingRhs (v) );
285  std::vector<Real> subSysCaRhs ( computeLocalSubSysCaRhs (v) );
286 
287  rhs[0] = computeLocalPotentialRhs (v );
288 
289  if (rhs.size() > M_numberOfEquations)
290  {
291  std::copy ( gatingRhs.begin(), gatingRhs.end(), rhs.begin() + 1 );
292  std::copy ( subSysCaRhs.begin(), subSysCaRhs.end(), rhs.begin() + 1 + gatingRhs.size() );
293  }
294  else
295  {
296  std::copy ( gatingRhs.begin(), gatingRhs.end(), rhs.begin() + 1 );
297  std::copy ( subSysCaRhs.begin(), subSysCaRhs.end() - 1, rhs.begin() + 1 + gatingRhs.size() );
298  }
299 
300 }
301 
302 void IonicFox::computeNonGatingRhs ( const std::vector<Real>& v, std::vector<Real>& rhs )
303 {
304  // std::vector<Real> subSysCaRhs ( computeLocalSubSysCaRhs (v) );
305  // std::copy ( subSysCaRhs.begin(), subSysCaRhs.end() - 1, rhs.begin() + M_numberOfGatingVariables );
306  std::vector<Real> gatingRhs ( computeLocalGatingRhs (v) );
307  std::vector<Real> subSysCaRhs ( computeLocalSubSysCaRhs (v) );
308 
309  std::copy ( gatingRhs.begin(), gatingRhs.end(), rhs.begin() );
310  std::copy ( subSysCaRhs.begin(), subSysCaRhs.end() - 1, rhs.begin() + gatingRhs.size() );
311 }
312 
313 Real IonicFox::computeLocalPotentialRhs ( const std::vector<Real>& v )
314 {
315  std::vector<Real> fastNa (fastINa (v) );
316  std::vector<Real> rapidK (rapidIK (v) );
317  std::vector<Real> slowK (slowIK (v) );
318  std::vector<Real> transOutK (transOutIK (v) );
319  std::vector<Real> subSysCaRHS (computeLocalSubSysCaRhs (v) );
320  Real INa = fastNa[0];
321  Real IKr = rapidK[0];
322  Real IKs = slowK[0];
323  Real It0 = transOutK[0];
324  Real ICa = subSysCaRHS[5];
325 
326  return - ( INa + timeIIK1 (v) + IKr + IKs + It0 + plaIKp (v) + pumpINaK (v)
327  + exINaCa (v) + backINab (v) + backICab (v) + pumpIpCa (v) + ICa ) ;
328 }
329 
330 std::vector<Real> IonicFox::computeLocalGatingRhs ( const std::vector<Real>& v )
331 {
332  Real m ( v[1] );
333  Real h ( v[2] );
334  Real j ( v[3] );
335 
336  Real xr ( v[4] );
337  Real xks ( v[5] );
338  Real xt0 ( v[6] );
339  Real yt0 ( v[7] );
340 
341  std::vector<Real> gatingRhs (7);
342  std::vector<Real> paramNa ( fastINa (v) );
343  std::vector<Real> paramKr ( rapidIK (v) );
344  std::vector<Real> paramKs ( slowIK (v) );
345  std::vector<Real> paramKt0 ( transOutIK (v) );
346 
347  gatingRhs[0] = paramNa[1] * ( 1 - m ) - paramNa[2] * m;
348  gatingRhs[1] = paramNa[3] * ( 1 - h ) - paramNa[4] * h;
349  gatingRhs[2] = paramNa[5] * ( 1 - j ) - paramNa[6] * j;
350 
351  gatingRhs[3] = ( paramKr[1] - xr ) / paramKr[2];
352  gatingRhs[4] = ( paramKs[1] - xks ) / paramKs[2];
353  gatingRhs[5] = paramKt0[1] * ( 1 - xt0 ) - paramKt0[2] * xt0;
354  gatingRhs[6] = paramKt0[3] * ( 1 - yt0 ) - paramKt0[4] * yt0;
355 
356  return gatingRhs;
357 }
358 
359 void IonicFox::computeGatingVariablesWithRushLarsen ( std::vector<Real>& v, const Real dt )
360 {
361  Real m ( v[1] );
362  Real h ( v[2] );
363  Real j ( v[3] );
364 
365  Real xr ( v[4] );
366  Real xks ( v[5] );
367  Real xt0 ( v[6] );
368  Real yt0 ( v[7] );
369 
370  std::vector<Real> paramNa ( fastINa (v) );
371  std::vector<Real> paramKr ( rapidIK (v) );
372  std::vector<Real> paramKs ( slowIK (v) );
373  std::vector<Real> paramKt0 ( transOutIK (v) );
374 
375  Real TAU_M = 1.0 / ( paramNa[1] + paramNa[2] );
376  Real M_INF = paramNa[1] * TAU_M;
377  Real TAU_H = 1.0 / ( paramNa[3] + paramNa[4] );
378  Real H_INF = paramNa[3] * TAU_M;
379  Real TAU_J = 1.0 / ( paramNa[5] + paramNa[6] );
380  Real J_INF = paramNa[5] * TAU_M;
381 
382  Real TAU_Xr = paramKr[2];
383  Real Xr_INF = paramKr[1];
384  Real TAU_Xks = paramKs[2];
385  Real Xks_INF = paramKs[1];
386 
387  Real TAU_xto = 1.0 / ( paramKt0[1] + paramKt0[2] );
388  Real xto_INF = paramKt0[1] * TAU_M;
389  Real TAU_yto = 1.0 / ( paramKt0[3] + paramKt0[4] );
390  Real yto_INF = paramKt0[3] * TAU_M;
391 
392  v[1] = M_INF - ( M_INF - m ) * std::exp (- dt / TAU_M );
393  v[2] = H_INF - ( H_INF - h ) * std::exp (- dt / TAU_H );
394  v[3] = J_INF - ( J_INF - j ) * std::exp (- dt / TAU_J );
395  v[4] = Xr_INF - ( Xr_INF - xr ) * std::exp (- dt / TAU_Xr );
396  v[5] = Xks_INF - ( Xks_INF - xks ) * std::exp (- dt / TAU_Xks );
397  v[6] = xto_INF - ( xto_INF - xt0 ) * std::exp (- dt / TAU_xto );
398  v[7] = yto_INF - ( yto_INF - yt0 ) * std::exp (- dt / TAU_yto );
399 }
400 
401 
402 //! Ca2+ Subsystem
403 
404 std::vector<Real> IonicFox::computeLocalSubSysCaRhs ( const std::vector<Real>& v )
405 {
406  std::vector<Real> subSysCaRHS (6);
407 
408  Real V ( v[0] );
409  Real cCaIn ( v[11] );
410  Real cCaSR ( v[12] );
411 
412  Real gating_d ( v[9] );
413  Real gating_f ( v[8] );
414  Real gating_f_Ca ( v[10] );
415 
416  // Internal Parameters
417 
418  Real gamma = 1.0 / ( 1.0 + pow ( ( 2000.0 / cCaSR ), 3) );
419 
420  Real jRel = M_Prel * gating_d * gating_f * gating_f_Ca * ( gamma * cCaSR - cCaIn ) / ( 1.0 + 1.65 * exp (V / 20.0) );
421  Real jLeak = M_Pleak * ( cCaSR - cCaIn );
422  Real jUp = M_Vup / ( 1.0 + pow ( ( M_KmUp / cCaIn ) , 2.0 ) );
423 
424  Real betaIn = 1 / ( 1 + M_CmdnTot * M_KmCmdn / ( ( M_KmCmdn + cCaIn ) * ( M_KmCmdn + cCaIn ) ) );
425  Real betaSR = 1 / ( 1 + M_CsqnTot * M_KmCsqn / ( ( M_KmCsqn + cCaSR ) * ( M_KmCsqn + cCaSR ) ) );
426 
427  Real sd = 1.0 / ( 1.0 + exp ( - ( V + 10.0 ) / 6.24 ) );
428  Real td = 1.0 / ( (0.25 * exp (-0.01 * V) / (1.0 + exp (-0.07 * V) ) ) + (0.07 * exp (-0.05 * (V + 40.0) ) / (1.0 + exp (0.05 * (V + 40.0) ) ) ) );
429 
430  Real sfca = 1.0 / (1.0 + pow ( (cCaIn / M_KmfCa), 3.0) );
431  Real tfca = 30.0;
432 
433  Real sf = 1.0 / (1 + exp (V + 12.5) / 5.0);
434  Real tf = 30.0 + (200.0 / (1 + exp ( (V + 20.0) / 9.5) ) );
435 
436  Real I_Ca_full = (M_PCa / M_Cm) * (4.0 * ( V * M_F * M_F ) / ( M_R * M_T ) ) *
437  ( (cCaIn * exp (2.0 * ( V * M_F ) / ( M_R * M_T ) ) - 0.341 * M_CaO) /
438  (exp (2.0 * ( V * M_F ) / ( M_R * M_T ) ) - 1.0) );
439  Real I_Cal = gating_d * gating_f * gating_f_Ca * I_Ca_full;
440  Real I_Ca_K = (M_PCaK / M_Cm) * (gating_d * gating_f * gating_f_Ca / (1.0 + (I_Ca_full / M_ICaHalf) ) ) *
441  (1000.0 * ( V * M_F * M_F ) / ( M_R * M_T ) ) * ( (M_KIn * exp ( ( V * M_F )
442  / ( M_R * M_T ) ) - M_K0) / (exp ( ( V * M_F ) / ( M_R * M_T ) ) - 1.0) );
443 
444  Real I_Ca = I_Cal + I_Ca_K;
445 
446  // RHS of the Ca2+ subsystem
447 
448  subSysCaRHS[4] = betaSR * ( jUp - jLeak - jRel) * ( M_VMyo / M_Vsr ); // dCaIn
449  subSysCaRHS[3] = betaIn * ( jRel + jLeak - jUp - (M_ACap * M_Cm / (2 * M_F * M_VMyo) ) * (I_Cal + backICab (v) + M_IpCa - 2 * exINaCa (v) ) ); //dCaSR
450 
451  subSysCaRHS[0] = ( sf - gating_f ) / tf;
452  subSysCaRHS[1] = ( sd - gating_d ) / td;
453  subSysCaRHS[2] = ( sfca - gating_f_Ca ) / tfca;
454 
455  subSysCaRHS[5] = I_Ca;
456 
457  return subSysCaRHS;
458 }
459 
460 //! Ionic Currents (Luo and Rudy)
461 
462 // Fast Na+ Current INa
463 std::vector<Real> IonicFox::fastINa ( const std::vector<Real>& v )
464 {
465  std::vector<Real> fastNa (7);
466 
467  Real V ( v[0] );
468  Real m ( v[1] );
469  Real h ( v[2] );
470  Real j ( v[3] );
471 
472  Real potNa = ( M_R * M_T / M_F ) * log ( M_NaO / M_NaIn );
473 
474  fastNa[0] = M_GNa * m * m * m * h * j * ( V - potNa );
475 
476  Real alpha_m = 0.32 * ( V + 47.13 ) / ( 1.0 - exp ( - 0.1 * ( V + 47.13 ) ) );
477  Real beta_m = 0.08 * exp (- V / 11.0 );
478 
479  Real alpha_h = 0.135 * exp ( - ( V + 80.0 ) / 6.8 );
480  Real beta_h = 7.5 / ( 1.0 + exp ( - 0.1 * ( V + 11.0 ) ) );
481  Real alpha_j = ( 0.175 * exp ( - ( V + 100.0 ) / 23.0 ) ) / ( 1.0 + exp ( 0.15 * ( V + 79.0 ) ) );
482  Real beta_j = 0.3 / ( 1.0 + exp ( -0.1 * ( V + 32.0 ) ) );
483 
484  fastNa[1] = alpha_m;
485  fastNa[2] = beta_m;
486  fastNa[3] = alpha_h;
487  fastNa[4] = beta_h;
488  fastNa[5] = alpha_j;
489  fastNa[6] = beta_j;
490 
491  return fastNa;
492 }
493 
494 // rapid K+ current
495 std::vector<Real> IonicFox::rapidIK ( const std::vector<Real>& v )
496 {
497  std::vector<Real> rapidK (3);
498 
499  Real V ( v[0] );
500  Real xr ( v[4] );
501 
502  Real potK1 = ( M_R * M_T / M_F ) * std::log ( M_K0 / M_KIn );
503  Real RV = 1.0 / ( 1.0 + 2.5 * std::exp ( 0.1 * ( V + 28.0 ) ) );
504 
505  Real sxr = 1.0 / ( 1.0 + std::exp ( -2.182 - 0.1819 * V ) );
506  Real txr = 43.0 + ( 1.0 / ( std::exp ( -5.495 + 0.1691 * V ) + std::exp ( -7.677 - 0.0128 * V ) ) );
507 
508  rapidK[0] = M_GKr * xr * RV * std::sqrt ( M_K0 / 4.0 ) * ( V - potK1 );
509 
510  rapidK[1] = sxr;
511  rapidK[2] = txr;
512 
513  return rapidK;
514 }
515 
516 // slow K+ current
517 std::vector<Real> IonicFox::slowIK ( const std::vector<Real>& v )
518 {
519 
520 
521  Real V ( v[0] );
522  Real xs ( v[5] );
523 
524  Real potKs = ( M_R * M_T / M_F ) * std::log ( ( M_K0 + 0.01833 * M_NaO ) / ( M_KIn + 0.01833 * M_NaIn ) );
525 
526  Real sxs = 1.0 / ( 1.0 + std::exp ( - ( V - 16.0 ) / 13.6 ) );
527  Real txs = 1.0 / ( ( 0.0000719 * ( V - 10.0 ) / ( 1.0 - std::exp ( -0.148 * ( V - 10.0 ) ) ) )
528  + ( 0.000131 * ( V - 10.0 ) / ( std::exp ( 0.0687 * (V - 10.0) ) - 1.0) ) );
529 
530  std::vector<Real> slow (3, 0.0);
531  slow[0] = M_GKs * xs * xs * ( V - potKs );
532  slow[1] = sxs;
533  slow[2] = txs;
534 
535  return slow;
536 }
537 
538 // transient outward K+ current
539 std::vector<Real> IonicFox::transOutIK ( const std::vector<Real>& v )
540 {
541  std::vector<Real> transOutK (5);
542 
543  Real V ( v[0] );
544  Real xto ( v[6] );
545  Real yto ( v[7] );
546 
547  Real potK1 = ( M_R * M_T / M_F ) * log ( M_K0 / M_KIn );
548 
549  Real axto = 0.04516 * exp ( 0.03577 * V );
550  Real bxto = 0.0989 * exp ( -0.06237 * V );
551  Real ayto = ( 0.005415 * exp ( - ( V + 33.5 ) / 5.0 ) ) / ( 1.0 + 0.051335 * exp ( - ( V + 33.5 ) / 5.0 ) );
552  Real byto = ( 0.005415 * exp ( ( V + 33.5 ) / 5.0 ) ) / ( 1.0 + 0.051335 * exp ( ( V + 33.5 ) / 5.0 ) );
553 
554  transOutK[0] = M_Gt0 * xto * yto * ( V - potK1 );
555 
556  transOutK[1] = axto;
557  transOutK[2] = bxto;
558  transOutK[3] = ayto;
559  transOutK[4] = byto;
560 
561  return transOutK;
562 }
563 
564 // Time-independent K+ Current IK1
565 Real IonicFox::timeIIK1 ( const std::vector<Real>& v )
566 {
567  Real V ( v[0] );
568 
569  Real potK1 = ( M_R * M_T / M_F ) * log ( M_K0 / M_KIn );
570  Real maxCondK1 = M_GK1;
571  Real sK1 = 1.0 / ( 2.0 + exp ( ( 1.62 / ( M_R * M_T / M_F ) ) * ( V - potK1 ) ) );
572  Real constK1inf = M_K0 / ( M_K0 + M_KmK1 );
573 
574  return maxCondK1 * sK1 * constK1inf * ( V - potK1 );
575 }
576 
577 // Plateau K+ current IKp
578 Real IonicFox::plaIKp ( const std::vector<Real>& v )
579 {
580  Real V ( v[0] );
581 
582  Real potKp = ( M_R * M_T / M_F ) * log ( M_K0 / M_KIn );
583  Real constKp = 1 / ( 1 + exp ( ( 7.488 - V ) / 5.98 ) );
584 
585  return M_GKp * constKp * ( V - potKp );
586 }
587 
588 // Na+/Ca2+ exchanger current INaCa
589 Real IonicFox::exINaCa ( const std::vector<Real>& v )
590 {
591 
592  Real V ( v[0] );
593  Real cCaIn ( v[11] );
594 
595  return M_kNaCa * ( 1.0 / ( pow (M_KmNa, 3) + pow (M_NaO, 3) ) ) * ( 1.0 / ( M_KmCa + M_CaO ) )
596  * (1.0 / ( 1.0 + M_kSat * exp ( ( M_eta - 1.0 ) * ( V * M_F ) / ( M_R * M_T ) ) ) )
597  * ( exp ( M_eta * ( V * M_F ) / ( M_R * M_T ) ) * pow (M_NaIn, 3) * M_CaO -
598  exp ( ( M_eta - 1.0 ) * ( V * M_F ) / ( M_R * M_T ) ) * pow (M_NaO, 3) * cCaIn );
599 }
600 
601 // Na+/K+ pump INaK
602 Real IonicFox::pumpINaK ( const std::vector<Real>& v )
603 {
604  Real V ( v[0] );
605 
606  Real sigma = ( 1.0 / 7.0 ) * ( exp ( M_NaO / 67.3 ) - 1.0 );
607  Real fNak = 1.0 / ( 1.0 + 0.1245 * exp ( -0.1 * ( V * M_F ) / ( M_R * M_T ) ) +
608  0.0365 * sigma * exp ( - ( V * M_F ) / ( M_R * M_T ) ) );
609 
610  return M_INaK * fNak * ( 1.0 / ( 1.0 + pow ( M_KmNai / M_NaIn, 1.5) ) ) * ( M_K0 / ( M_K0 + M_KmK0 ) );
611 }
612 
613 // Sarcolemmal Ca2+ pump current IpCa
614 Real IonicFox::pumpIpCa ( const std::vector<Real>& v )
615 {
616  return M_IpCa * v[11] / ( M_KmPCa + v[11] );
617 }
618 
619 // Ca2+ background current ICab OK
620 Real IonicFox::backICab ( const std::vector<Real>& v )
621 {
622  Real V ( v[0] );
623  Real cCaIn ( v[11] );
624 
625  Real potCaN = ( ( M_R * M_T ) / ( 2 * M_F ) ) * log ( M_CaO / cCaIn );
626 
627  return M_GCab * ( V - potCaN );
628 }
629 
630 // Na+ background current INab
631 Real IonicFox::backINab ( const std::vector<Real>& v)
632 {
633  Real V ( v[0] );
634 
635  Real potNaN = M_R * M_T / M_F * log ( M_NaO / M_NaIn );
636 
637  return M_GNab * ( V - potNaN );
638 }
639 
640 void IonicFox::showMe()
641 {
642  std::cout << "\n\n\t\tIonicFox Informations\n\n";
643  std::cout << "number of unkowns: " << this->Size() << std::endl;
644 
645  // std::cout << "\n\t\tList of model parameters:\n\n";
646  // std::cout << "areaCap: " << this->areaCap() << std::endl;
647  // std::cout << "volMyo: " << this->volMyo() << std::endl;
648  // std::cout << "volJSR: " << this->volJSR() << std::endl;
649  // std::cout << "volNSR: " << this->volNSR() << std::endl;
650  // std::cout << "volSS: " << this->volSS() << std::endl;
651  // std::cout << "concNa0: " << this->concNa0() << std::endl;
652  // std::cout << "concCa0: " << this->concCa0() << std::endl;
653  // std::cout << "lTrpnTot: " << this->lTrpnTot() << std::endl;
654  // std::cout << "hTrpnTot: " << this->hTrpnTot() << std::endl;
655  // std::cout << "kpHtrpn: " << this->kpHtrpn() << std::endl;
656  // std::cout << "knHtrpn: " << this->knHtrpn() << std::endl;
657  // std::cout << "kpLtrpn: " << this->kpLtrpn() << std::endl;
658  // std::cout << "knLtrpn: " << this->knLtrpn() << std::endl;
659  // std::cout << "cmdnTot: " << this->cmdnTot() << std::endl;
660  // std::cout << "csqnTot: " << this->csqnTot() << std::endl;
661  // std::cout << "constmCmdn: " << this->constmCmdn() << std::endl;
662  // std::cout << "constmCsqn: " << this->constmCsqn() << std::endl;
663  // std::cout << "capMem: " << this->capMem() << std::endl;
664  // std::cout << "farad: " << this->farad() << std::endl;
665  // std::cout << "temp: " << this->temp() << std::endl;
666  // std::cout << "gasConst: " << this->gasConst() << std::endl;
667  // std::cout << "maxCondNa: " << this->maxCondNa() << std::endl;
668  // std::cout << "maxCondKp: " << this->maxCondKp() << std::endl;
669  // std::cout << "permNaK: " << this->permNaK() << std::endl;
670  // std::cout << "KNaCa: " << this->kNaCa() << std::endl;
671  // std::cout << "constmNa: " << this->constmNa() << std::endl;
672  // std::cout << "constmCa: " << this->constmCa() << std::endl;
673  // std::cout << "kSat: " << this->kSat() << std::endl;
674  // std::cout << "eta: " << this->eta() << std::endl;
675  // std::cout << "courNaK: " << this->courNaK() << std::endl;
676  // std::cout << "constmNai: " << this->constmNai() << std::endl;
677  // std::cout << "constmK0: " << this->constmK0() << std::endl;
678  // std::cout << "permNsCa: " << this->permNsCa() << std::endl;
679  // std::cout << "constmNsCa: " << this->constmNsCa() << std::endl;
680  // std::cout << "courpCa: " << this->courpCa() << std::endl;
681  // std::cout << "constmpCa: " << this->constmCa() << std::endl;
682  // std::cout << "maxCondCab: " << this->maxCondCab() << std::endl;
683  // std::cout << "maxCondNab: " << this->maxCondNab() << std::endl;
684  // std::cout << "maxRyRPerm: " << this->maxRyRPerm() << std::endl;
685  // std::cout << "leakRateConst: " << this->leakRateConst() << std::endl;
686  // std::cout << "pumpRateATPase: " << this->pumpRateATPase() << std::endl;
687  // std::cout << "constmUp: " << this->constmUp() << std::endl;
688  // std::cout << "timeConstNsrJsr: " << this->timeConstNsrJsr() << std::endl;
689  // std::cout << "timeConstSubMyo: " << this->timeConstSubMyo() << std::endl;
690  // std::cout << "kAPlus: " << this->kAPlus() << std::endl;
691  // std::cout << "kANeg: " << this->kANeg() << std::endl;
692  // std::cout << "kBPlus: " << this->kBPlus() << std::endl;
693  // std::cout << "kBNeg: " << this->kBNeg() << std::endl;
694  // std::cout << "kCPlus: " << this->kCPlus() << std::endl;
695  // std::cout << "kCNeg: " << this->kCNeg() << std::endl;
696  // std::cout << "coopParamN: " << this->coopParamN() << std::endl;
697  // std::cout << "coopParamM: " << this->coopParamM() << std::endl;
698  // std::cout << "intoOpenSt: " << this->intoOpenSt() << std::endl;
699  // std::cout << "outOpenSt: " << this->outOpenSt() << std::endl;
700  // std::cout << "intoOpenStca: " << this->intoOpenStCa() << std::endl;
701  // std::cout << "outOpentSt2: " << this->outOpenSt2() << std::endl;
702  // std::cout << "modeTParamA: " << this->modeTParamA() << std::endl;
703  // std::cout << "modeTParamB: " << this->modeTParamB() << std::endl;
704  // std::cout << "modeTParamO: " << this->modeTParamO() << std::endl;
705  // std::cout << "permCa: " << this->permCa() << std::endl;
706  // std::cout << "permK: " << this->permK() << std::endl;
707  // std::cout << "courCaHalf: " << this->courCaHalf() << std::endl;
708 
709 
710  std::cout << "\n\t\t End of IonicFox Informations\n\n\n";
711 }
712 
713 
714 }
LifeV::IonicFox::timeIIK1
Real timeIIK1(const std::vector< Real > &v)
Definition: IonicFox.cpp:565
LifeV::IonicFox::exINaCa
Real exINaCa(const std::vector< Real > &v)
Definition: IonicFox.cpp:589
LifeV::IonicFox::showMe
void showMe()
Display information about the model.
Definition: IonicFox.cpp:640
LifeV::IonicFox::M_ACap
Real M_ACap
Cell Geometry Parameters (4)
Definition: IonicFox.hpp:550
LifeV::IonicFox::super
ElectroIonicModel super
Definition: IonicFox.hpp:67
LifeV::IonicFox::M_PCa
Real M_PCa
L-type Ca2+ Channel Parameters (6)
Definition: IonicFox.hpp:603
LifeV::IonicFox::computeLocalGatingRhs
std::vector< Real > computeLocalGatingRhs(const std::vector< Real > &v)
Definition: IonicFox.cpp:330
LifeV::IonicFox::IonicFox
IonicFox()
Constructor.
Definition: IonicFox.cpp:51
LifeV::IonicFox::M_NaO
Real M_NaO
Standard Ionic Concentrations (5)
Definition: IonicFox.hpp:557
LifeV::IonicFox::slowIK
std::vector< Real > slowIK(const std::vector< Real > &v)
Definition: IonicFox.cpp:517
LifeV::IonicFox::computeRhs
void computeRhs(const std::vector< Real > &v, std::vector< Real > &rhs)
This methods contains the actual evaluation of the rhs of all state variablesin the model (0D version...
Definition: IonicFox.cpp:282
LifeV::IonicFox::M_Cm
Real M_Cm
Membrane Current Parameters (24)
Definition: IonicFox.hpp:572
LifeV::IonicFox::rapidIK
std::vector< Real > rapidIK(const std::vector< Real > &v)
Definition: IonicFox.cpp:495
ETCurrentFE::updateInverseJacobian
void updateInverseJacobian(const UInt &iQuadPt)
Definition: ETCurrentFE.cpp:405
LifeV::IonicFox::M_CmdnTot
Real M_CmdnTot
Buffering Parameters (4)
Definition: IonicFox.hpp:565
LifeV::IonicFox::pumpIpCa
Real pumpIpCa(const std::vector< Real > &v)
Definition: IonicFox.cpp:614
LifeV::IonicFox::computeGatingVariablesWithRushLarsen
void computeGatingVariablesWithRushLarsen(std::vector< Real > &v, const Real dt)
This methods contains the actual evaluation of the rhs of the voltage equation only (0D version) ...
Definition: IonicFox.cpp:359
LifeV::IonicFox::IonicFox
IonicFox(Teuchos::ParameterList &parameterList)
Definition: IonicFox.cpp:112
LifeV::IonicFox::M_KmUp
Real M_KmUp
SR Parameters (1)
Definition: IonicFox.hpp:599
LifeV::IonicFox::plaIKp
Real plaIKp(const std::vector< Real > &v)
Definition: IonicFox.cpp:578
LifeV::IonicFox::computeLocalSubSysCaRhs
std::vector< Real > computeLocalSubSysCaRhs(const std::vector< Real > &v)
Ca2+ Subsystem.
Definition: IonicFox.cpp:404
LifeV::IonicFox
IonicModel - This class implements an ionic model.
Definition: IonicFox.hpp:61
LifeV::IonicFox::operator=
IonicFox & operator=(const IonicFox &model)
Operator.
Definition: IonicFox.cpp:214
LifeV::Real
double Real
Generic real data.
Definition: LifeV.hpp:175
LifeV::IonicFox::backICab
Real backICab(const std::vector< Real > &v)
Definition: IonicFox.cpp:620
LifeV::IonicFox::backINab
Real backINab(const std::vector< Real > &v)
Definition: IonicFox.cpp:631
LifeV::IonicFox::transOutIK
std::vector< Real > transOutIK(const std::vector< Real > &v)
Definition: IonicFox.cpp:539
LifeV::IonicFox::computeGatingRhs
void computeGatingRhs(const std::vector< Real > &v, std::vector< Real > &rhs)
Methods.
Definition: IonicFox.cpp:272
LifeV::IonicFox::pumpINaK
Real pumpINaK(const std::vector< Real > &v)
Definition: IonicFox.cpp:602
LifeV::IonicFox::computeLocalPotentialRhs
Real computeLocalPotentialRhs(const std::vector< Real > &v)
This methods contains the actual evaluation of the rhs of the voltage equation only (0D version) ...
Definition: IonicFox.cpp:313
LifeV::IonicFox::computeNonGatingRhs
void computeNonGatingRhs(const std::vector< Real > &v, std::vector< Real > &rhs)
In the case this method is improperly used, it should use this default implementation.
Definition: IonicFox.cpp:302
LifeV::IonicFox::IonicFox
IonicFox(const IonicFox &model)
Definition: IonicFox.cpp:160
LifeV::IonicFox::fastINa
std::vector< Real > fastINa(const std::vector< Real > &v)
Ionic Currents (Luo and Rudy)
Definition: IonicFox.cpp:463