LifeV
IonicHodgkinHuxley.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 Hodgkin and Huxley
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/IonicHodgkinHuxley.hpp>
40 
41 
42 namespace LifeV
43 {
44 
45 IonicHodgkinHuxley::IonicHodgkinHuxley() :
46  super ( 4, 3 ),
47  M_gNa (120.0),
48  M_gK (36.0),
49  M_gL (0.3),
50  M_vNa (115.0),
51  M_vK (-12.0),
52  M_vL (10.6)
53 {
54  M_restingConditions.at (0) = 0.0;
55  M_restingConditions.at (1) = 0.052932485257250;
56  M_restingConditions.at (2) = 0.317676914060697;
57  M_restingConditions.at (3) = 0.596120753508460;
58 
59 }
60 
61 IonicHodgkinHuxley::IonicHodgkinHuxley ( Teuchos::ParameterList& parameterList ) :
62  super ( 4, 3 )
63 {
64  M_gNa = parameterList.get ("gNa", 120.0 );
65  M_gK = parameterList.get ("gK", 36.0 );
66  M_gL = parameterList.get ("gL", 0.3 );
67  M_vNa = parameterList.get ("vNa", 115.0 );
68  M_vK = parameterList.get ("vK", -12.0 );
69  M_vL = parameterList.get ("vL", 10.6 );
70 
71  M_restingConditions.at (0) = 0.0;
72  M_restingConditions.at (1) = 0.052932485257250;
73  M_restingConditions.at (2) = 0.317676914060697;
74  M_restingConditions.at (3) = 0.596120753508460;
75 }
76 
77 IonicHodgkinHuxley::IonicHodgkinHuxley ( const IonicHodgkinHuxley& model )
78 {
79 
80  M_gNa = model.M_gNa;
81  M_gK = model.M_gK;
82  M_gL = model.M_gL;
83  M_vNa = model.M_vNa;
84  M_vK = model.M_vK;
85  M_vL = model.M_vL;
86 
87 
88  M_numberOfEquations = model.M_numberOfEquations;
89  M_restingConditions = model.M_restingConditions;
90  M_numberOfGatingVariables = model.M_numberOfGatingVariables;
91 
92 }
93 
94 // ===================================================
95 //! Operator
96 // ===================================================
97 IonicHodgkinHuxley& IonicHodgkinHuxley::operator= ( const IonicHodgkinHuxley& model )
98 {
99  M_gNa = model.M_gNa;
100  M_gK = model.M_gK;
101  M_gL = model.M_gL;
102  M_vNa = model.M_vNa;
103  M_vK = model.M_vK;
104  M_vL = model.M_vL;
105 
106  M_numberOfEquations = model.M_numberOfEquations;
107  M_restingConditions = model.M_restingConditions;
108  M_numberOfGatingVariables = model.M_numberOfGatingVariables;
109 
110 
111  return *this;
112 }
113 
114 
115 // ===================================================
116 //! Methods
117 // ===================================================
118 void IonicHodgkinHuxley::computeGatingRhs ( const std::vector<Real>& v,
119  std::vector<Real>& rhs )
120 {
121 
122  Real V = v[0];
123  Real M = v[1];
124  Real N = v[2];
125  Real H = v[3];
126 
127  Real alpham = 0.1 * (25. - V) / (std::exp ( (25. - V) / 10.) - 1.);
128  Real betam = 4.*std::exp (-V / 18.0);
129  Real alphah = 0.07 * std::exp (-V / 20.);
130  Real betah = 1.0 / (std::exp ( (30. - V) / 10.) + 1.);
131  Real alphan = 0.01 * (10. - V) / (std::exp ( (10. - V) / 10.) - 1.);
132  Real betan = 0.125 * std::exp (-V / 80.0);
133 
134  rhs[0] = alpham * (1 - M) - betam * M;
135  rhs[1] = alphan * (1 - N) - betan * N;
136  rhs[2] = alphah * (1 - H) - betah * H;
137 }
138 
139 void IonicHodgkinHuxley::computeRhs ( const std::vector<Real>& v,
140  std::vector<Real>& rhs )
141 {
142  std::vector<Real> tmpRhs (this->Size() - 1, 0.0);
143  computeGatingRhs (v, tmpRhs);
144  rhs[0] = computeLocalPotentialRhs (v);
145  rhs[1] = tmpRhs[0];
146  rhs[2] = tmpRhs[1];
147  rhs[3] = tmpRhs[2];
148 }
149 
150 
151 Real IonicHodgkinHuxley::computeLocalPotentialRhs ( const std::vector<Real>& v )
152 {
153  Real dPotential (0.0);
154 
155  Real V = v[0];
156  Real M = v[1];
157  Real N = v[2];
158  Real H = v[3];
159 
160 
161  dPotential = -M_gK * N * N * N * N * (V - M_vK) - M_gNa * M * M * M * H * (V - M_vNa) - M_gL * (V - M_vL); //+M_appliedCurrent;
162 
163  return dPotential;
164 }
165 
166 void IonicHodgkinHuxley::computeGatingVariablesWithRushLarsen ( std::vector<Real>& v, const Real dt )
167 {
168  Real V = v[0];
169  Real M = v[1];
170  Real N = v[2];
171  Real H = v[3];
172 
173  Real alpham = 0.1 * (25. - V) / (std::exp ( (25. - V) / 10.) - 1.);
174  Real betam = 4.*std::exp (-V / 18.0);
175  Real alphah = 0.07 * std::exp (-V / 20.);
176  Real betah = 1.0 / (std::exp ( (30. - V) / 10.) + 1.);
177  Real alphan = 0.01 * (10. - V) / (std::exp ( (10. - V) / 10.) - 1.);
178  Real betan = 0.125 * std::exp (-V / 80.0);
179 
180  Real taum = alpham + betam;
181  Real taun = alphan + betan;
182  Real tauh = alphah + betah;
183 
184  Real mInf = alpham / (taum);
185  Real nInf = alphan / (taun);
186  Real hInf = alphah / (tauh);
187 
188  v[1] = mInf + (M - mInf) * std::exp (-dt * taum);
189  v[2] = nInf + (N - nInf) * std::exp (-dt * taun);
190  v[3] = hInf + (H - hInf) * std::exp (-dt * tauh);
191 
192 }
193 
194 void IonicHodgkinHuxley::showMe()
195 {
196 
197  std::cout << "\n\tHi, I'm the Hodgkin Huxley model for neurons. This is the first model implemented by Simone Palamara!!!\n\t See you soon\n\n";
198 }
199 
200 
201 }
ETCurrentFE::updateInverseJacobian
void updateInverseJacobian(const UInt &iQuadPt)
Definition: ETCurrentFE.cpp:405
LifeV::Real
double Real
Generic real data.
Definition: LifeV.hpp:175