LifeV
IonicLuoRudyI.hpp
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 IonicLuoRudyI
29  @brief Ionic model Luo-Rudy I
30 
31  model as in
32  Luo, Ching-hsing, and Yoram Rudy.
33  "A model of the ventricular cardiac action potential.
34  Depolarization, repolarization, and their interaction."
35  Circulation research 68.6 (1991): 1501-1526.
36 
37  @date 01-2013
38  @author Simone Rossi <simone.rossi@epfl.ch>
39 
40  @contributors
41  @mantainer Simone Rossi <simone.rossi@epfl.ch>
42  @last update 01-2013
43  */
44 
45 
46 #ifndef _IONICLUORUDYI_H_
47 #define _IONICLUORUDYI_H_
48 
49 #include <lifev/electrophysiology/solver/IonicModels/ElectroIonicModel.hpp>
50 
51 #include <Teuchos_RCP.hpp>
52 #include <Teuchos_ParameterList.hpp>
53 #include "Teuchos_XMLParameterListHelpers.hpp"
54 
55 namespace LifeV
56 {
57 //! IonicModel - This class implements an ionic model.
58 
59 class IonicLuoRudyI : public virtual ElectroIonicModel
60 {
61 
62 public:
63  //! @name Type definitions
64  //@{
65  typedef ElectroIonicModel super;
66  typedef std::shared_ptr<VectorEpetra> vectorPtr_Type;
67  typedef std::shared_ptr<VectorElemental> elvecPtr_Type;
68  typedef RegionMesh<LinearTetra> mesh_Type;
69  //@}
70 
71 
72 
73  //! @name Constructors & Destructor
74  //@{
75 
76  //! Constructor
77  IonicLuoRudyI();
78 
79  /*!
80  * @param Epetra communicator
81  * @param list of parameters in an xml file
82  */
83  IonicLuoRudyI ( Teuchos::ParameterList& parameterList );
84 
85  /*!
86  * @param IonicLuoRudyI object
87  */
88  IonicLuoRudyI ( const IonicLuoRudyI& model );
89  //! Destructor
90  virtual ~IonicLuoRudyI() {}
91 
92  //@}
93 
94  //! @name Overloads
95  //@{
96 
97  IonicLuoRudyI& operator= ( const IonicLuoRudyI& model );
98 
99  //@}
100 
101  //! @name Setters and getters
102  //@{
103 
104  //parameters getters and setters
105  inline const Real ENa() const
106  {
107  return M_ENa;
108  }
109  inline const Real gNa() const
110  {
111  return M_gNa;
112  }
113 
114  // //Slow Inward Current
115  inline const Real gsi() const
116  {
117  return M_gsi;
118  }
119 
120  // //Time dependent potassium current
121  inline const Real K0() const
122  {
123  return M_K0;
124  }
125  inline const Real gK() const
126  {
127  return M_gK;
128  }
129  inline const Real EK() const
130  {
131  return M_EK;
132  }
133 
134  // //Time independent potassium current
135  inline const Real gK1() const
136  {
137  return M_gK1;
138  }
139  inline const Real EK1() const
140  {
141  return M_EK1;
142  }
143 
144  // //Plateau Current
145  inline const Real EKp() const
146  {
147  return M_EKp;
148  }
149  inline const Real gKp() const
150  {
151  return M_gKp;
152  }
153 
154  // //Background Current
155  inline const Real gb() const
156  {
157  return M_gb;
158  }
159 
160 
161  //set methods
162  inline void setENa ( const Real p )
163  {
164  M_ENa = p;
165  }
166  inline void setGNa ( const Real p )
167  {
168  M_gNa = p;
169  }
170 
171  // //Slow Inward Current
172  inline void setGsi ( const Real p )
173  {
174  M_gsi = p;
175  }
176 
177  // //Time dependent potassium current
178  inline void setK0 ( const Real p )
179  {
180  M_K0 = p;
181  M_gK = computeGK (M_K0);
182  M_gK1 = computeGK1 (M_K0);
183  }
184  inline void setGK ( const Real p )
185  {
186  M_gK = p;
187  }
188  inline void setEK ( const Real p )
189  {
190  M_EK = p;
191  }
192 
193  // //Time independent potassium current
194  inline void setGK1 ( const Real p )
195  {
196  M_gK1 = p;
197  }
198  inline void setEK1 ( const Real p )
199  {
200  M_EK1 = p;
201  }
202 
203  // //Plateau Current
204  inline void setEKp ( const Real p )
205  {
206  M_EKp = p;
207  }
208  inline void setGKp ( const Real p )
209  {
210  M_gKp = p;
211  }
212 
213  // //Background Current
214  inline void setGb ( const Real p )
215  {
216  M_gb = p;
217  }
218 
219  //inline const short int& Size() const { return M_numberOfEquations; }
220  //@}
221 
222  //! @name Methods
223  //@{
224 
225  //Fast inward current
226  //===================
227  //gating variable m
228  inline Real am ( Real V )
229  {
230  return ( 0.32 * ( V + 47.13 ) / ( 1.0 - std::exp ( - 0.1 * ( V + 47.13 ) ) ) );
231  }
232  inline Real bm ( Real V )
233  {
234  return ( 0.08 * std::exp ( - V / 11.0 ) );
235  }
236  inline Real tm (Real V)
237  {
238  return ( 1.0 / ( am (V) + bm (V) ) );
239  }
240  inline Real minf (Real V)
241  {
242  return ( am (V) * tm (V) );
243  }
244  inline Real dm (Real V, Real m)
245  {
246  return ( am (V) * (1.0 - m ) - bm (V) * m );
247  }
248 
249  //gating variable m
250  inline Real ah ( Real V )
251  {
252  if ( V < - 40.0 )
253  {
254  return ( 0.135 * std::exp ( - ( V + 80 ) / 6.8 ) );
255  }
256  else
257  {
258  return 0.0;
259  }
260  }
261  inline Real bh ( Real V )
262  {
263  if ( V < - 40.0 )
264  {
265  return ( ( 3.56 * std::exp ( 0.079 * V ) + 3.1 * 1e5 * std::exp ( 0.35 * V) ) );
266  }
267  else
268  {
269  return ( 1.0 / ( 0.13 * ( 1 + std::exp ( - ( V + 10.66 ) / 11.1 ) ) ) );
270  }
271  }
272  inline Real th (Real V)
273  {
274  return ( 1.0 / ( ah (V) + bh (V) ) );
275  }
276  inline Real hinf (Real V)
277  {
278  return ( ah (V) * th (V) );
279  }
280  inline Real dh (Real V, Real h)
281  {
282  return ( ah (V) * (1.0 - h ) - bh (V) * h );
283  }
284 
285  //gating variable j
286  inline Real aj ( Real V )
287  {
288  if ( V < - 40.0 )
289  return ( ( -1.2714 * 1e5 * std::exp ( 0.2444 * V ) - 3.474 * 1e-5 * std::exp ( - 0.04391 * V ) )
290  * ( V + 37.78 ) / ( 1 + std::exp ( 0.311 * ( V + 79.23 ) ) ) );
291  else
292  {
293  return 0.0;
294  }
295  }
296  inline Real bj ( Real V )
297  {
298  if ( V < - 40.0 )
299  {
300  return ( 0.1212 * std::exp ( - 0.01052 * V ) / ( 1 + std::exp ( -0.1378 * ( V + 40.14 ) ) ) );
301  }
302  else
303  {
304  return ( 0.3 * std::exp ( -2.535 * 1e-7 * V ) / ( 1 + std::exp ( -0.1 * ( V + 32 ) ) ) );
305  }
306  }
307  inline Real tj (Real V)
308  {
309  return ( 1.0 / ( aj (V) + bj (V) ) );
310  }
311  inline Real jinf (Real V)
312  {
313  return ( aj (V) * tj (V) );
314  }
315  inline Real dj (Real V, Real j)
316  {
317  return ( aj (V) * (1.0 - j ) - bj (V) * j );
318  }
319 
320  //Fast Inward Current
321  inline Real INa (Real V, Real m, Real h, Real j)
322  {
323  return ( M_gNa * m * m * m * h * j * ( V - M_ENa ) );
324  }
325 
326  // // slow inward current %%
327  // //=========================
328  inline Real Esi ( Real Ca )
329  {
330  return ( 7.7 - 13.0287 * std::log (Ca) );
331  }
332  //gating variable d
333  inline Real ad ( Real V )
334  {
335  return ( 0.095 * std::exp ( - 0.01 * ( V - 5 ) ) / ( 1 + std::exp ( - 0.072 * ( V - 5 ) ) ) );
336  }
337  inline Real bd ( Real V )
338  {
339  return ( 0.07 * std::exp ( - 0.017 * ( V + 44) ) / ( 1 + std::exp ( 0.05 * ( V + 44) ) ) );
340  }
341  inline Real td (Real V)
342  {
343  return ( 1.0 / ( ad (V) + bd (V) ) );
344  }
345  inline Real dinf (Real V)
346  {
347  return ( ad (V) * td (V) );
348  }
349  inline Real dd (Real V, Real d)
350  {
351  return ( ad (V) * (1.0 - d ) - bd (V) * d );
352  }
353 
354  //gating variable f
355  inline Real af ( Real V )
356  {
357  return ( 0.012 * std::exp ( - 0.008 * ( V + 28 ) ) / ( 1 + std::exp ( 0.15 * ( V + 28 ) ) ) );
358  }
359  inline Real bf ( Real V )
360  {
361  return ( 0.0065 * std::exp ( - 0.02 * ( V + 30 ) ) / ( 1 + std::exp ( - 0.2 * ( V + 30 ) ) ) );
362  }
363  inline Real tf (Real V)
364  {
365  return ( 1.0 / ( af (V) + bf (V) ) );
366  }
367  inline Real finf (Real V)
368  {
369  return ( af (V) * tf (V) );
370  }
371  inline Real df (Real V, Real f)
372  {
373  return ( af (V) * (1.0 - f ) - bf (V) * f );
374  }
375 
376  //Slow inward current
377  inline Real Isi (Real V, Real d, Real f, Real Ca)
378  {
379  return ( M_gsi * d * f * ( V - Esi (Ca) ) );
380  }
381 
382  //Calcium concentration
383  inline Real dCa (Real V, Real d, Real f, Real Ca)
384  {
385  return ( - 1e-4 * Isi (V, d, f, Ca) + 0.07 * ( 1e-4 - Ca ) );
386  }
387 
388  //Time dependent potassium current %%
389  //===================================
390  inline Real computeGK (Real K0)
391  {
392  return ( 0.282 * std::sqrt ( K0 / 5.4) );
393  }
394 
395  //gating variable X
396  inline Real aX ( Real V )
397  {
398  return ( 0.0005 * std::exp ( 0.083 * ( V + 50 ) ) / ( 1 + std::exp ( 0.057 * ( V + 50 ) ) ) );
399  }
400  inline Real bX ( Real V )
401  {
402  return ( 0.0013 * std::exp ( - 0.06 * ( V + 20 ) ) / ( 1 + std::exp ( - 0.04 * ( V + 20 ) ) ) );
403  }
404  inline Real tX (Real V)
405  {
406  return ( 1.0 / ( aX (V) + bX (V) ) );
407  }
408  inline Real Xinf (Real V)
409  {
410  return ( aX (V) * tX (V) );
411  }
412  inline Real dX (Real V, Real X)
413  {
414  return ( aX (V) * (1.0 - X ) - bX (V) * X );
415  }
416 
417  //gating Xi
418  inline Real Xi (Real V)
419  {
420  if ( V > -100 )
421  {
422  return ( 2.837 * ( std::exp ( 0.04 * ( V + 77 ) ) - 1 ) / ( ( V + 77 ) * std::exp ( 0.04 * ( V + 35 ) ) ) );
423  }
424  else
425  {
426  return 1.0;
427  }
428  }
429  //Time dependent potassium current
430  inline Real IK (Real V, Real X)
431  {
432  return M_gK * X * Xi (V) * ( V - M_EK );
433  }
434 
435  //Time independent potassium current
436  //===================================
437  inline Real computeGK1 (Real K0)
438  {
439  return ( 0.6047 * std::sqrt ( K0 / 5.4) );
440  }
441 
442  //gating variable K1
443  inline Real aK1 ( Real V )
444  {
445  return ( 1.02 / ( 1 + std::exp ( 0.2385 * ( V - M_EK1 - 59.2915 ) ) ) );
446  }
447  inline Real bK1 ( Real V )
448  {
449  return ( ( 0.49124 * std::exp ( 0.08032 * ( V - M_EK1 + 5.4760 ) )
450  + exp ( 0.06175 * ( V - M_EK1 - 594.31 ) ) ) / ( 1 + std::exp ( - 0.5143 * ( V - M_EK1 + 4.753 ) ) ) );
451  }
452  inline Real tK1 (Real V)
453  {
454  return ( 1.0 / ( aK1 (V) + bK1 (V) ) );
455  }
456  inline Real K1inf (Real V)
457  {
458  return ( aK1 (V) * tK1 (V) );
459  }
460 
461  //Time independent potassium current
462  inline Real IK1 (Real V)
463  {
464  return M_gK1 * K1inf (V) * ( V - M_EK1 );
465  }
466 
467  //Plateau Current
468  //====================
469  inline Real Kp ( Real V )
470  {
471  return ( 1.0 / ( 1 + std::exp ( ( 7.488 - V ) / 5.98 ) ) );
472  }
473  inline Real IKp (Real V)
474  {
475  return M_gKp * Kp (V) * ( V - M_EKp );
476  }
477 
478  //Background Current
479  //====================
480  inline Real Ib (Real V)
481  {
482  return M_gb * ( V + 59.87 );
483  }
484 
485  //Total current
486  //=============
487  inline Real Itot (Real V, Real m, Real h, Real j, Real d, Real f, Real X, Real Ca)
488  {
489  return ( INa (V, m, h, j) + Isi (V, d, f, Ca) + IK (V, X) + IK1 (V) + IKp (V) + Ib (V) );
490  }
491  //Compute the rhs on a single node or for the 0D case
492  void computeGatingRhs ( const std::vector<Real>& v, std::vector<Real>& rhs);
493 
494  void computeNonGatingRhs ( const std::vector<Real>& v, std::vector<Real>& rhs );
495 
496 
497  void computeRhs ( const std::vector<Real>& v, std::vector<Real>& rhs);
498 
499  void computeGatingVariablesWithRushLarsen ( std::vector<Real>& v, const Real dt );
500 
501  // compute the rhs with state variable interpolation
502  Real computeLocalPotentialRhs ( const std::vector<Real>& v );
503 
504  //! Display information about the model
505  void showMe();
506 
507  //! Solves the ionic model
508  //virtual void solveXbModel( const vector_Type& Calcium,
509  // const Real timeStep )=0;
510  //@}
511 
512 private:
513  //! Model Parameters
514 
515  //! Chemical kinetics parameters
516  //Fast Sodium Current
517  Real M_ENa;
518  Real M_gNa;
519 
520  //Slow Inward Current
521  Real M_gsi;
522 
523  //Time dependent potassium current
524  Real M_K0;
525  Real M_gK;
526  Real M_EK;
527 
528  //Time independent potassium current
529  Real M_gK1;
530  Real M_EK1;
531 
532  //Plateau Current
533  Real M_EKp;
534  Real M_gKp;
535  //Background Current
536  Real M_gb;
537  //! Xb states == equivalent to the number of equations
538  //short int M_numberOfEquations;
539 
540  //@}
541 
542 }; // class IonicLuoRudyI
543 
544 
545 
546 }
547 
548 #endif
ETCurrentFE::updateInverseJacobian
void updateInverseJacobian(const UInt &iQuadPt)
Definition: ETCurrentFE.cpp:405
LifeV::Real
double Real
Generic real data.
Definition: LifeV.hpp:175