LifeV
IonicMinimalModel.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 Bueno-Orovio et al.
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/IonicMinimalModel.hpp>
40 #include <lifev/eta/fem/ETFESpace.hpp>
41 #include <lifev/eta/expression/Integrate.hpp>
42 #include <boost/typeof/typeof.hpp>
43 
44 namespace LifeV
45 {
46 
47 // ===================================================
48 //! Constructors
49 // ===================================================
50 IonicMinimalModel::IonicMinimalModel() :
51  super ( 4, 2 ),
52  M_uo ( 0. ),
53  M_uu ( 1.58 ),
54  M_tetav ( 0.3 ),
55  M_tetaw ( 0.015 ),
56  M_tetavm ( 0.015 ),
57  M_tetao ( 0.006 ),
58  M_tauv1 ( 60.0 ),
59  M_tauv2 ( 1150.0),
60  M_tauvp ( 1.4506),
61  M_tauw1 ( 70.0 ),
62  M_tauw2 ( 20.0 ),
63  M_kw ( 65.0 ),
64  M_uw ( 0.03 ),
65  M_tauwp ( 280.0 ),
66  M_taufi ( 0.11 ),
67  M_tauo1 ( 6.0 ),
68  M_tauo2 ( 6.0 ),
69  M_tauso1 ( 43.0 ),
70  M_tauso2 ( 0.2 ),
71  M_kso ( 2.0 ),
72  M_uso ( 0.65 ),
73  M_taus1 ( 2.7342),
74  M_taus2 ( 3.0 ),
75  M_ks ( 2.0994),
76  M_us ( 0.9087),
77  M_tausi ( 2.8723),
78  M_tauwinf ( 0.07 ),
79  M_winfstar ( 0.94 )
80 {
81  M_restingConditions.at (0) = 0.0;
82  M_restingConditions.at (1) = 1.0;
83  M_restingConditions.at (2) = 1.0;
84  M_restingConditions.at (3) = 0.021553043080281;
85 
86 }
87 
88 IonicMinimalModel::IonicMinimalModel ( Teuchos::ParameterList& parameterList ) :
89  super ( 4 )
90 {
91  M_uo = parameterList.get ("uo", 1000.0 );
92  M_uu = parameterList.get ("uu", 1.61 );
93  M_tetav = parameterList.get ("tetav", 0.30 );
94  M_tetaw = parameterList.get ("tetaw", 0.130 );
95  M_tetavm = parameterList.get ("tetavm", 0.10 );
96  M_tetao = parameterList.get ("tetao", 0.005 );
97  M_tauv1 = parameterList.get ("tauv1", 80.0 );
98  M_tauv2 = parameterList.get ("tauv2", 1.4506 );
99  M_tauvp = parameterList.get ("tauvp", 1.4506 );
100  M_tauw1 = parameterList.get ("tauw1", 70.0 );
101  M_tauw2 = parameterList.get ("tauw2", 8.0 );
102  M_kw = parameterList.get ("kw", 200.0 );
103  M_uw = parameterList.get ("uw", 0.016 );
104  M_tauwp = parameterList.get ("tauwp", 280.0 );
105  M_taufi = parameterList.get ("taufi", 0.078 );
106  M_tauo1 = parameterList.get ("tauo1", 410.0 );
107  M_tauo2 = parameterList.get ("tauo2", 7.0 );
108  M_tauso1 = parameterList.get ("tauso1", 91.0 );
109  M_tauso2 = parameterList.get ("tauso2", 0.8 );
110  M_kso = parameterList.get ("kso", 2.1 );
111  M_uso = parameterList.get ("uso", 0.6 );
112  M_taus1 = parameterList.get ("taus1", 2.7342 );
113  M_taus2 = parameterList.get ("taus2", 4.0 );
114  M_ks = parameterList.get ("ks", 2.0994 );
115  M_us = parameterList.get ("us", 0.9087 );
116  M_tausi = parameterList.get ("tausi", 3.3849 );
117  M_tauwinf = parameterList.get ("tauwinf", 0.01 );
118  M_winfstar = parameterList.get ("winfstar", 0.5 );
119 
120  M_restingConditions.at (0) = 0.0;
121  M_restingConditions.at (1) = 1.0;
122  M_restingConditions.at (2) = 1.0;
123  M_restingConditions.at (3) = 0.021553043080281;
124 }
125 
126 IonicMinimalModel::IonicMinimalModel ( const IonicMinimalModel& model )
127 {
128 
129  M_uo = model.M_uo;
130  M_uu = model.M_uu;
131  M_tetav = model.M_tetav;
132  M_tetaw = model.M_tetaw;
133  M_tetavm = model.M_tetavm;
134  M_tetao = model.M_tetao;
135  M_tauv1 = model.M_tauv1;
136  M_tauv2 = model.M_tauv2;
137  M_tauvp = model.M_tauvp;
138  M_tauw1 = model.M_tauw1;
139  M_tauw2 = model.M_tauw2;
140  M_kw = model.M_kw;
141  M_uw = model.M_uw;
142  M_tauwp = model.M_tauwp;
143  M_taufi = model.M_taufi;
144  M_tauo1 = model.M_tauo1;
145  M_tauo2 = model.M_tauo2;
146  M_tauso1 = model.M_tauso1;
147  M_tauso2 = model.M_tauso2;
148  M_kso = model.M_kso;
149  M_uso = model.M_uso;
150  M_taus1 = model.M_taus1;
151  M_taus2 = model.M_taus2;
152  M_ks = model.M_ks;
153  M_us = model.M_us;
154  M_tausi = model.M_tausi;
155  M_tauwinf = model.M_tauwinf;
156  M_winfstar = model.M_winfstar;
157 
158  M_numberOfEquations = model.M_numberOfEquations;
159  M_restingConditions = model.M_restingConditions;
160 
161 }
162 
163 // ===================================================
164 //! Operator
165 // ===================================================
166 IonicMinimalModel& IonicMinimalModel::operator= ( const IonicMinimalModel& model )
167 {
168  M_uo = model.M_uo;
169  M_uu = model.M_uu;
170  M_tetav = model.M_tetav;
171  M_tetaw = model.M_tetaw;
172  M_tetavm = model.M_tetavm;
173  M_tetao = model.M_tetao;
174  M_tauv1 = model.M_tauv1;
175  M_tauv2 = model.M_tauv2;
176  M_tauvp = model.M_tauvp;
177  M_tauw1 = model.M_tauw1;
178  M_tauw2 = model.M_tauw2;
179  M_kw = model.M_kw;
180  M_uw = model.M_uw;
181  M_tauwp = model.M_tauwp;
182  M_taufi = model.M_taufi;
183  M_tauo1 = model.M_tauo1;
184  M_tauo2 = model.M_tauo2;
185  M_tauso1 = model.M_tauso1;
186  M_tauso2 = model.M_tauso2;
187  M_kso = model.M_kso;
188  M_uso = model.M_uso;
189  M_taus1 = model.M_taus1;
190  M_taus2 = model.M_taus2;
191  M_ks = model.M_ks;
192  M_us = model.M_us;
193  M_tausi = model.M_tausi;
194  M_tauwinf = model.M_tauwinf;
195  M_winfstar = model.M_winfstar;
196 
197  M_numberOfEquations = model.M_numberOfEquations;
198  M_restingConditions = model.M_restingConditions;
199 
200 
201  return *this;
202 }
203 
204 
205 // ===================================================
206 //! Methods
207 // ===================================================
208 void IonicMinimalModel::computeGatingRhs ( const std::vector<Real>& v,
209  std::vector<Real>& rhs )
210 {
211 
212  Real U = v[0];
213  Real V = v[1];
214  Real W = v[2];
215  Real S = v[3];
216 
217  Real tauvm = ( 1.0 - Heaviside ( U - M_tetavm ) ) * M_tauv1 + Heaviside ( U - M_tetavm ) * M_tauv2;
218  Real tauwm = M_tauw1 + ( M_tauw2 - M_tauw1 ) * ( 1.0 + std::tanh ( M_kw * ( U - M_uw ) ) ) / 2.0;
219  Real taus = ( 1.0 - Heaviside ( U - M_tetaw ) ) * M_taus1 + Heaviside ( U - M_tetaw ) * M_taus2;
220 
221  Real vinf = Heaviside ( M_tetavm - U );
222  Real winf = ( 1.0 - Heaviside ( U - M_tetao ) ) * ( 1.0 - U / M_tauwinf ) + Heaviside ( U - M_tetao ) * M_winfstar;
223 
224  rhs[0] = ( 1.0 - Heaviside ( U - M_tetav ) ) * ( vinf - V ) / tauvm - Heaviside ( U - M_tetav ) * V / M_tauvp;
225  rhs[1] = ( 1.0 - Heaviside ( U - M_tetaw ) ) * ( winf - W ) / tauwm - Heaviside ( U - M_tetaw ) * W / M_tauwp;
226  rhs[2] = ( ( 1.0 + std::tanh ( M_ks * ( U - M_us ) ) ) / 2.0 - S ) / taus;
227 
228 }
229 
230 void IonicMinimalModel::computeRhs ( const std::vector<Real>& v,
231  std::vector<Real>& rhs )
232 {
233 
234  Real U = v[0];
235  Real V = v[1];
236  Real W = v[2];
237  Real S = v[3];
238 
239  Real tauvm = ( 1.0 - Heaviside ( U - M_tetavm ) ) * M_tauv1 + Heaviside ( U - M_tetavm ) * M_tauv2;
240  Real tauwm = M_tauw1 + ( M_tauw2 - M_tauw1 ) * ( 1.0 + std::tanh ( M_kw * ( U - M_uw ) ) ) / 2.0;
241  Real tauso = M_tauso1 + ( M_tauso2 - M_tauso1 ) * ( 1.0 + std::tanh ( M_kso * ( U - M_uso ) ) ) / 2.0;
242  Real taus = ( 1.0 - Heaviside ( U - M_tetaw ) ) * M_taus1 + Heaviside ( U - M_tetaw ) * M_taus2;
243  Real tauo = ( 1.0 - Heaviside ( U - M_tetao ) ) * M_tauo1 + Heaviside ( U - M_tetao ) * M_tauo2;
244 
245  Real vinf = Heaviside ( M_tetavm - U );
246  Real winf = ( 1.0 - Heaviside ( U - M_tetao ) ) * ( 1.0 - U / M_tauwinf ) + Heaviside ( U - M_tetao ) * M_winfstar;
247 
248  Real Jfi = - V * Heaviside ( U - M_tetav ) * ( U - M_tetav ) * ( M_uu - U ) / M_taufi;
249  Real Jso = ( U - M_uo ) * ( 1.0 - Heaviside ( U - M_tetaw ) ) / tauo + Heaviside ( U - M_tetaw ) / tauso;
250  Real Jsi = - Heaviside ( U - M_tetaw ) * W * S / M_tausi;
251 
252  rhs[0] = - ( Jfi + Jso + Jsi );
253  rhs[1] = ( 1.0 - Heaviside ( U - M_tetav ) ) * ( vinf - V ) / tauvm - Heaviside ( U - M_tetav ) * V / M_tauvp;
254  rhs[2] = ( 1.0 - Heaviside ( U - M_tetaw ) ) * ( winf - W ) / tauwm - Heaviside ( U - M_tetaw ) * W / M_tauwp;
255  rhs[3] = ( ( 1.0 + std::tanh ( M_ks * ( U - M_us ) ) ) / 2.0 - S ) / taus;
256 
257 }
258 
259 
260 void IonicMinimalModel::computeGatingVariablesWithRushLarsen ( std::vector<Real>& v, const Real dt )
261 {
262  std::cout << "\n\nRush Larsen method, for minimal model not implemented!!!\n";
263  std::cout << "\n\nI will use Forward Euler!!!\n";
264  std::vector<Real> rhs (3, 0.0);
265  computeGatingRhs ( v, rhs );
266  v[1] += dt * rhs[0];
267  v[2] += dt * rhs[1];
268  v[3] += dt * rhs[2];
269 }
270 
271 
272 Real IonicMinimalModel::computeLocalPotentialRhs ( const std::vector<Real>& v )
273 {
274  Real dPotential (0.0);
275 
276  Real U = v[0];
277  Real V = v[1];
278  Real W = v[2];
279  Real S = v[3];
280 
281  Real tauso = M_tauso1 + ( M_tauso2 - M_tauso1 ) * ( 1.0 + std::tanh ( M_kso * ( U - M_uso ) ) ) / 2.0;
282  Real tauo = ( 1.0 - Heaviside ( U - M_tetao ) ) * M_tauo1 + Heaviside ( U - M_tetao ) * M_tauo2;
283 
284  Real Jfi = - V * Heaviside ( U - M_tetav ) * ( U - M_tetav ) * ( M_uu - U ) / M_taufi;
285  Real Jso = ( U - M_uo ) * ( 1.0 - Heaviside ( U - M_tetaw ) ) / tauo + Heaviside ( U - M_tetaw ) / tauso;
286  Real Jsi = - Heaviside ( U - M_tetaw ) * W * S / M_tausi;
287 
288  dPotential = - ( Jfi + Jso + Jsi );
289  // std::cout.precision(15);
290  // std::cout << "\nValue: " << Jso;
291  return dPotential;
292 }
293 
294 
295 void IonicMinimalModel::showMe()
296 {
297 
298  std::cout << "\n\tHi, I'm the minimal model\n\t See you soon\n\n";
299 
300  std::cout << "\nuo " << M_uo ;
301  std::cout << "\nuu " << M_uu ;
302  std::cout << "\ntetav " << M_tetav ;
303  std::cout << "\ntetaw " << M_tetaw ;
304  std::cout << "\ntetavm " << M_tetavm ;
305  std::cout << "\ntetao " << M_tetao ;
306  std::cout << "\ntauv1 " << M_tauv1 ;
307  std::cout << "\ntauv2 " << M_tauv2 ;
308  std::cout << "\ntauvp " << M_tauvp ;
309  std::cout << "\ntauw1 " << M_tauw1 ;
310  std::cout << "\ntauw2 " << M_tauw2 ;
311  std::cout << "\nkw " << M_kw ;
312  std::cout << "\nuw " << M_uw ;
313  std::cout << "\ntauwp " << M_tauwp ;
314  std::cout << "\ntaufi " << M_taufi ;
315  std::cout << "\ntauo1 " << M_tauo1 ;
316  std::cout << "\ntauo2 " << M_tauo2 ;
317  std::cout << "\ntauso1 " << M_tauso1 ;
318  std::cout << "\ntauso2 " << M_tauso2 ;
319  std::cout << "\nkso " << M_kso ;
320  std::cout << "\nuso " << M_uso ;
321  std::cout << "\ntaus1 " << M_taus1 ;
322  std::cout << "\ntaus2 " << M_taus2 ;
323  std::cout << "\nks " << M_ks ;
324  std::cout << "\nus " << M_us ;
325  std::cout << "\ntausi " << M_tausi ;
326  std::cout << "\ntauwinf " << M_tauwinf ;
327  std::cout << "\nwinfstar " << M_winfstar << "\n" ;
328 }
329 
330 void IonicMinimalModel::computePotentialRhsSVI ( const std::vector<vectorPtr_Type>& v,
331  std::vector<vectorPtr_Type>& rhs,
332  FESpace<mesh_Type, MapEpetra>& uFESpace)
333 {
334  typedef ETFESpace<mesh_Type, MapEpetra, 3, 1> ETFESpace_Type;
335  typedef std::shared_ptr<ETFESpace_Type> ETFESpacePtr_Type;
336 
337  * (rhs[0]) *= 0.0;
338 
339  if (uFESpace.mapPtr() -> commPtr() -> MyPID() == 0)
340  {
341  std::cout << "\nMinimal Model: Assembling SVI using ETA!\n";
342  }
343 
344  ETFESpacePtr_Type spaceScalar (
345  new ETFESpace_Type (uFESpace.mesh(), &feTetraP1, uFESpace.mapPtr() -> commPtr() ) );
346 
347  {
348  using namespace ExpressionAssembly;
349 
350  std::shared_ptr<MMTanhFunctor> tanh (new MMTanhFunctor);
351  std::shared_ptr<MMHFunctor> H (new MMHFunctor);
352  std::shared_ptr<MMSV> sv (new MMSV);
353 
354  BOOST_AUTO_TPL (U, value (spaceScalar, * (v[0] ) ) );
355  BOOST_AUTO_TPL (V, value (spaceScalar, * (v[1] ) ) );
356  BOOST_AUTO_TPL (W, value (spaceScalar, * (v[2] ) ) );
357  BOOST_AUTO_TPL (S, value (spaceScalar, * (v[3] ) ) );
358  BOOST_AUTO_TPL (Iapp, value (spaceScalar, *M_appliedCurrentPtr ) );
359 
360 
361  BOOST_AUTO_TPL (tauso, M_tauso1 + ( M_tauso2 - M_tauso1 ) * ( 1.0 + eval (tanh, M_kso * ( U - M_uso ) ) ) / 2.0);
362  BOOST_AUTO_TPL (tauo, ( 1.0 - eval (H, U - M_tetao ) ) * M_tauo1 + eval (H, U - M_tetao ) * M_tauo2);
363  BOOST_AUTO_TPL (Jfi, value (-1.0) * V * eval (H, U - M_tetav ) * ( U - M_tetav ) * ( M_uu - U ) / M_taufi);
364  BOOST_AUTO_TPL (Jso, ( U - M_uo ) * ( 1.0 - eval (H, U - M_tetaw ) ) / tauo + eval (H, U - M_tetaw ) / tauso);
365  BOOST_AUTO_TPL (Jsi, value (-1.0) * eval (H, U - M_tetaw ) * W * S / M_tausi);
366 
367  integrate ( elements ( uFESpace.mesh() ),
368  uFESpace.qr(),
369  spaceScalar,
370  ( Iapp - ( Jfi + Jso + Jsi ) ) * phi_i ) >> rhs.at (0);
371  }
372 
373  rhs.at (0) -> globalAssemble();
374 }
375 
376 }