Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
FentonKarma.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File FentonKarma.cpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) n6 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 // License for the specific language governing rights and limitations under
13 // Permission is hereby granted, free of charge, to any person obtaining a
14 // copy of this software and associated documentation files (the "Software"),
15 // to deal in the Software without restriction, including without limitation
16 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 // and/or sell copies of the Software, and to permit persons to whom the
18 // Software is furnished to do so, subject to the following conditions:
19 //
20 // The above copyright notice and this permission notice shall be included
21 // in all copies or substantial portions of the Software.
22 //
23 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 // DEALINGS IN THE SOFTWARE.
30 //
31 // Description: Fenton-Karma cell model.
32 //
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 #include <iostream>
36 #include <string>
37 
38 #include <LibUtilities/BasicUtils/VmathArray.hpp>
39 #include <CardiacEPSolver/CellModels/FentonKarma.h>
40 
41 using namespace std;
42 
43 namespace Nektar
44 {
45  // Register cell model
46  std::string FentonKarma::className =
47  GetCellModelFactory().RegisterCreatorFunction(
48  "FentonKarma",
49  FentonKarma::create,
50  "Phenomenological Model.");
51 
52  // Register cell model variants
53  std::string FentonKarma::lookupIds[10] = {
54  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
55  "BR", FentonKarma::eBR),
56  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
57  "MBR", FentonKarma::eMBR),
58  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
59  "MLR1", FentonKarma::eMLR1),
60  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
61  "GP", FentonKarma::eGP),
62  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
63  "CF1", FentonKarma::eCF1),
64  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
65  "CF2a", FentonKarma::eCF2a),
66  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
67  "CF2b", FentonKarma::eCF2b),
68  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
69  "CF2c", FentonKarma::eCF2c),
70  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
71  "CF3a", FentonKarma::eCF3a),
72  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
73  "CF3b", FentonKarma::eCF3b)
74  };
75 
76  // Register default variant
77  std::string FentonKarma::def =
78  LibUtilities::SessionReader::RegisterDefaultSolverInfo(
79  "CellModelVariant","BR");
80 
81  /**
82  * Initialise Fenton-Karma model parameters.
83  * k1 is k in the original model.
84  * k2 is an additional parameter introduced in Cherry&Fenton 2004.
85  * u_r and u_fi are introduced by Cherry&Fenton 2004 and is the same as
86  * u_c in the original model.
87  */
88  FentonKarma::FentonKarma(
89  const LibUtilities::SessionReaderSharedPtr& pSession,
90  const MultiRegions::ExpListSharedPtr& pField)
91  : CellModel(pSession, pField)
92  {
93  model_variant = pSession->GetSolverInfoAsEnum<FentonKarma::Variants>(
94  "CellModelVariant");
95 
96  C_m = 1; // picoF
97  V_0 = -85;
98 
99  switch (model_variant) {
100  case eBR:
101  g_fi_max = 4;
102  tau_r = 33.33;
103  tau_si = 29;
104  tau_0 = 12.5;
105  tau_v_plus = 3.33;
106  tau_v1_minus = 1250;
107  tau_v2_minus = 19.6;
108  tau_w_plus = 870;
109  tau_w_minus = 41;
110  u_c = 0.13;
111  u_v = 0.04;
112  u_r = 0.13;
113  u_fi = 0.13;
114  u_csi = 0.85;
115  k1 = 10;
116  k2 = 0.0;
117  break;
118  case eMBR:
119  g_fi_max = 4;
120  tau_r = 50;
121  tau_si = 44.84;
122  tau_0 = 8.3;
123  tau_v_plus = 3.33;
124  tau_v1_minus = 1000;
125  tau_v2_minus = 19.2;
126  tau_w_plus = 667;
127  tau_w_minus = 11;
128  u_c = 0.13;
129  u_v = 0.04;
130  u_r = 0.13;
131  u_fi = 0.13;
132  u_csi = 0.85;
133  k1 = 10;
134  k2 = 0.0;
135  break;
136  case eMLR1:
137  g_fi_max = 5.8;
138  tau_r = 130;
139  tau_si = 127;
140  tau_0 = 12.5;
141  tau_v_plus = 10;
142  tau_v1_minus = 18.2;
143  tau_v2_minus = 18.2;
144  tau_w_plus = 1020;
145  tau_w_minus = 80;
146  u_c = 0.13;
147  u_v = 0.0;
148  u_r = 0.13;
149  u_fi = 0.13;
150  u_csi = 0.85;
151  k1 = 10;
152  k2 = 0.0;
153  break;
154  case eGP:
155  g_fi_max = 8.7;
156  tau_r = 25;
157  tau_si = 22.22;
158  tau_0 = 12.5;
159  tau_v_plus = 10;
160  tau_v1_minus = 333;
161  tau_v2_minus = 40;
162  tau_w_plus = 1000;
163  tau_w_minus = 65;
164  u_c = 0.13;
165  u_v = 0.025;
166  u_r = 0.13;
167  u_fi = 0.13;
168  u_csi = 0.85;
169  k1 = 10;
170  k2 = 0.0;
171  break;
172  case eCF1:
173  g_fi_max = 6.6666;
174  tau_r = 12.5;
175  tau_si = 10;
176  tau_0 = 1.5;
177  tau_v_plus = 10;
178  tau_v1_minus = 350;
179  tau_v2_minus = 80;
180  tau_w_plus = 562;
181  tau_w_minus = 48.5;
182  u_c = 0.25;
183  u_v = 0.001;
184  u_r = 0.25;
185  u_fi = 0.15;
186  u_csi = 0.2;
187  k1 = 15;
188  k2 = 0;
189  break;
190  case eCF2a:
191  g_fi_max = 6.6666;
192  tau_r = 31;
193  tau_si = 26.5;
194  tau_0 = 1.5;
195  tau_v_plus = 10;
196  tau_v1_minus = 20;
197  tau_v2_minus = 20;
198  tau_w_plus = 800;
199  tau_w_minus = 45;
200  u_c = 0.25;
201  u_v = 0.05;
202  u_r = 0.6;
203  u_fi = 0.11;
204  u_csi = 0.7;
205  k1 = 10;
206  k2 = 1;
207  break;
208  case eCF2b:
209  g_fi_max = 6.6666;
210  tau_r = 31;
211  tau_si = 26.5;
212  tau_0 = 1.5;
213  tau_v_plus = 10;
214  tau_v1_minus = 100;
215  tau_v2_minus = 20;
216  tau_w_plus = 800;
217  tau_w_minus = 45;
218  u_c = 0.25;
219  u_v = 0.05;
220  u_r = 0.6;
221  u_fi = 0.11;
222  u_csi = 0.7;
223  k1 = 10;
224  k2 = 1;
225  break;
226  case eCF2c:
227  g_fi_max = 6.6666;
228  tau_r = 31;
229  tau_si = 26.5;
230  tau_0 = 1.5;
231  tau_v_plus = 10;
232  tau_v1_minus = 150;
233  tau_v2_minus = 20;
234  tau_w_plus = 800;
235  tau_w_minus = 45;
236  u_c = 0.25;
237  u_v = 0.05;
238  u_r = 0.6;
239  u_fi = 0.11;
240  u_csi = 0.7;
241  k1 = 10;
242  k2 = 1;
243  break;
244  case eCF3a:
245  g_fi_max = 13.3333;
246  tau_r = 38;
247  tau_si = 127;
248  tau_0 = 8.3;
249  tau_v_plus = 3.33;
250  tau_v1_minus = 45;
251  tau_v2_minus = 300;
252  tau_w_plus = 600;
253  tau_w_minus = 40;
254  tau_y_plus = 1000;
255  tau_y_minus = 230;
256  u_c = 0.25;
257  u_v = 0.5;
258  u_r = 0.25;
259  u_fi = 0.25;
260  u_csi = 0.7;
261  k1 = 60;
262  k2 = 0;
263  break;
264  case eCF3b:
265  g_fi_max = 13.3333;
266  tau_r = 38;
267  tau_si = 127;
268  tau_0 = 8.3;
269  tau_v_plus = 3.33;
270  tau_v1_minus = 20;
271  tau_v2_minus = 300;
272  tau_w_plus = 600;
273  tau_w_minus = 40;
274  tau_y_plus = 1000;
275  tau_y_minus = 230;
276  u_c = 0.25;
277  u_v = 0.5;
278  u_r = 0.25;
279  u_fi = 0.25;
280  u_csi = 0.7;
281  k1 = 60;
282  k2 = 0;
283  break;
284  }
285 
286  tau_d = C_m/g_fi_max;
287 
288  isCF3 = (model_variant == eCF3a || model_variant == eCF3b);
289 
290  // List gates and concentrations
291  m_gates.push_back(1);
292  m_gates.push_back(2);
293  m_nvar = 3;
294 
295  // Cherry-Fenton 2004 Model 3 has extra gating variable
296  if (isCF3)
297  {
298  m_gates.push_back(3);
299  m_nvar = 4;
300  }
301 
302  ASSERTL0(!isCF3, "Cherry-Fenton model 3 not implemented yet.");
303  }
304 
305 
306 
307  /**
308  *
309  */
310  FentonKarma::~FentonKarma()
311  {
312 
313  }
314 
315 
316 
317  void FentonKarma::v_Update(
318  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
319  Array<OneD, Array<OneD, NekDouble> >&outarray,
320  const NekDouble time )
321  {
322  ASSERTL0(inarray.get() != outarray.get(),
323  "Must have different arrays for input and output.");
324 
325  // Variables
326  // 0 u membrane potential
327  // 1 v v gate
328  // 2 w w gate
329  int n = m_nq;
330  int i = 0;
331 
332  // Declare pointers
333  const NekDouble *u = &inarray[0][0];
334  const NekDouble *v = &inarray[1][0];
335  const NekDouble *w = &inarray[2][0];
336  const NekDouble *y = isCF3 ? &inarray[3][0] : 0;
337  NekDouble *u_new = &outarray[0][0];
338  NekDouble *v_new = &outarray[1][0];
339  NekDouble *w_new = &outarray[2][0];
340  //NekDouble *y_new = isCF3 ? &outarray[3][0] : 0;
341  NekDouble *v_tau = &m_gates_tau[0][0];
342  NekDouble *w_tau = &m_gates_tau[1][0];
343  //NekDouble *y_tau = isCF3 ? &m_gates_tau[2][0] : 0;
344 
345  // Temporary variables
346  NekDouble J_fi, J_so, J_si, h1, h2, h3, alpha, beta;
347 
348  // Compute rates for each point in domain
349  for (i = 0; i < n; ++i)
350  {
351  // Heavyside functions
352  h1 = (*u < u_c) ? 0.0 : 1.0;
353  h2 = (*u < u_v) ? 0.0 : 1.0;
354  h3 = (*u < u_r) ? 0.0 : 1.0;
355 
356  // w-gate
357  alpha = (1-h1)/tau_w_minus;
358  beta = h1/tau_w_plus;
359  *w_tau = 1.0 / (alpha + beta);
360  *w_new = alpha * (*w_tau);
361 
362  // v-gate
363  alpha = (1-h1)/(h2*tau_v1_minus + (1-h2)*tau_v2_minus);
364  beta = h1/tau_v_plus;
365  *v_tau = 1.0 / (alpha + beta);
366  *v_new = alpha * (*v_tau);
367 
368  // y-gate
369  if (isCF3)
370  {
371  // TODO: implementation for y_tau and y_new
372  }
373 
374  // J_fi
375  J_fi = -(*v)*h1*(1 - *u)*(*u - u_fi)/tau_d;
376 
377  // J_so
378  // added extra (1-k2*v) term from Cherry&Fenton 2004
379  J_so = (*u)*(1-h3)*(1-k2*(*v))/tau_0 +
380  (isCF3 ? h3*(*u)*(*y)/tau_r : h3/tau_r);
381 
382  // J_si
383  J_si = -(*w)*(1 + tanh(k1*(*u - u_csi)))/(2.0*tau_si);
384 
385  // u
386  *u_new = -J_fi - J_so - J_si;
387 
388  ++u, ++v, ++w, ++u_new, ++v_new, ++w_new, ++v_tau, ++w_tau;
389  }
390  }
391 
392  void FentonKarma::v_GenerateSummary(SummaryList& s)
393  {
394  SolverUtils::AddSummaryItem(s, "Cell model", "Fenton Karma");
395  SolverUtils::AddSummaryItem(s, "Cell model var.", lookupIds[model_variant]);
396  }
397 
398 
399  void FentonKarma::v_SetInitialConditions()
400  {
401  Vmath::Fill(m_nq, 0.0, m_cellSol[0], 1);
402  Vmath::Fill(m_nq, 1.0, m_cellSol[1], 1);
403  Vmath::Fill(m_nq, 1.0, m_cellSol[2], 1);
404  if (isCF3)
405  {
406  Vmath::Fill(m_nq, 0.1, m_cellSol[2], 1);
407  }
408 
409  }
410 
411  std::string FentonKarma::v_GetCellVarName(unsigned int idx)
412  {
413  switch (idx)
414  {
415  case 0: return "u";
416  case 1: return "v";
417  case 2: return "w";
418  case 3: return "y";
419  default: return "unknown";
420  }
421  }
422 }
Generated on Sun Mar 15 2015 00:11:31 for Nektar++ by   doxygen 1.8.1.2