Nektar++
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 // Permission is hereby granted, free of charge, to any person obtaining a
13 // copy of this software and associated documentation files (the "Software"),
14 // to deal in the Software without restriction, including without limitation
15 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 // and/or sell copies of the Software, and to permit persons to whom the
17 // Software is furnished to do so, subject to the following conditions:
18 //
19 // The above copyright notice and this permission notice shall be included
20 // in all copies or substantial portions of the Software.
21 //
22 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 // DEALINGS IN THE SOFTWARE.
29 //
30 // Description: Fenton-Karma cell model.
31 //
32 ///////////////////////////////////////////////////////////////////////////////
33 
34 #include <iostream>
35 #include <string>
36 
37 #include <LibUtilities/BasicUtils/VmathArray.hpp>
38 #include <CardiacEPSolver/CellModels/FentonKarma.h>
39 
40 using namespace std;
41 
42 namespace Nektar
43 {
44  // Register cell model
45  std::string FentonKarma::className =
46  GetCellModelFactory().RegisterCreatorFunction(
47  "FentonKarma",
48  FentonKarma::create,
49  "Phenomenological Model.");
50 
51  // Register cell model variants
52  std::string FentonKarma::lookupIds[] = {
53  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
54  "BR", FentonKarma::eBR),
55  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
56  "MBR", FentonKarma::eMBR),
57  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
58  "MLR1", FentonKarma::eMLR1),
59  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
60  "GP", FentonKarma::eGP),
61  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
62  "CF1", FentonKarma::eCF1),
63  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
64  "CF2a", FentonKarma::eCF2a),
65  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
66  "CF2b", FentonKarma::eCF2b),
67  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
68  "CF2c", FentonKarma::eCF2c),
69  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
70  "CF3a", FentonKarma::eCF3a),
71  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
72  "CF3b", FentonKarma::eCF3b),
73  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
74  "FC2002Set1a", FentonKarma::eFC2002Set1a),
75  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
76  "FC2002Set1b", FentonKarma::eFC2002Set1b),
77  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
78  "FC2002Set1c", FentonKarma::eFC2002Set1c),
79  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
80  "FC2002Set1d", FentonKarma::eFC2002Set1d),
81  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
82  "FC2002Set1e", FentonKarma::eFC2002Set1e),
83  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
84  "FC2002Set2", FentonKarma::eFC2002Set2),
85  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
86  "FC2002Set4", FentonKarma::eFC2002Set4),
87  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
88  "FC2002Set5", FentonKarma::eFC2002Set5),
89  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
90  "FC2002Set6", FentonKarma::eFC2002Set6),
91  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
92  "FC2002Set7", FentonKarma::eFC2002Set7),
93  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
94  "FC2002Set8", FentonKarma::eFC2002Set8),
95  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
96  "FC2002Set9", FentonKarma::eFC2002Set9),
97  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
98  "Lawson", FentonKarma::eLawson),
99  LibUtilities::SessionReader::RegisterEnumValue("CellModelVariant",
100  "CAF", FentonKarma::eCAF)
101  };
102 
103  // Register default variant
104  std::string FentonKarma::def =
105  LibUtilities::SessionReader::RegisterDefaultSolverInfo(
106  "CellModelVariant","BR");
107 
108  /**
109  * Initialise Fenton-Karma model parameters.
110  * k1 is k in the original model.
111  * k2 is an additional parameter introduced in Cherry&Fenton 2004.
112  * u_r and u_fi are introduced by Cherry&Fenton 2004 and is the same as
113  * u_c in the original model.
114  */
115  FentonKarma::FentonKarma(
116  const LibUtilities::SessionReaderSharedPtr& pSession,
117  const MultiRegions::ExpListSharedPtr& pField)
118  : CellModel(pSession, pField)
119  {
120  model_variant = pSession->GetSolverInfoAsEnum<FentonKarma::Variants>(
121  "CellModelVariant");
122 
123  C_m = 1; // picoF
124  V_0 = -85;
125  V_fi = 15;
126 
127  switch (model_variant) {
128  case eBR: // (Fenton, Karma, Chaos 8(1), 1998)
129  g_fi_max = 4;
130  tau_r = 33.33;
131  tau_si = 29;
132  tau_0 = 12.5;
133  tau_v_plus = 3.33;
134  tau_v1_minus = 1250;
135  tau_v2_minus = 19.6;
136  tau_w_plus = 870;
137  tau_w_minus = 41;
138  u_c = 0.13;
139  u_v = 0.04;
140  u_r = 0.13;
141  u_fi = 0.13;
142  u_csi = 0.85;
143  k1 = 10;
144  k2 = 0.0;
145  break;
146  case eMBR: // (Fenton, Karma, Chaos 8(1), 1998)
147  g_fi_max = 4;
148  tau_r = 50;
149  tau_si = 44.84;
150  tau_0 = 8.3;
151  tau_v_plus = 3.33;
152  tau_v1_minus = 1000;
153  tau_v2_minus = 19.2;
154  tau_w_plus = 667;
155  tau_w_minus = 11;
156  u_c = 0.13;
157  u_v = 0.04;
158  u_r = 0.13;
159  u_fi = 0.13;
160  u_csi = 0.85;
161  k1 = 10;
162  k2 = 0.0;
163  break;
164  case eMLR1: // (Fenton, Karma, Chaos 8(1), 1998)
165  g_fi_max = 5.8;
166  tau_r = 130;
167  tau_si = 127;
168  tau_0 = 12.5;
169  tau_v_plus = 10;
170  tau_v1_minus = 18.2;
171  tau_v2_minus = 18.2;
172  tau_w_plus = 1020;
173  tau_w_minus = 80;
174  u_c = 0.13;
175  u_v = 0.0;
176  u_r = 0.13;
177  u_fi = 0.13;
178  u_csi = 0.85;
179  k1 = 10;
180  k2 = 0.0;
181  break;
182  case eGP: // (Fenton, Karma, Chaos 8(1), 1998)
183  g_fi_max = 8.7;
184  tau_r = 25;
185  tau_si = 22.22;
186  tau_0 = 12.5;
187  tau_v_plus = 10;
188  tau_v1_minus = 333;
189  tau_v2_minus = 40;
190  tau_w_plus = 1000;
191  tau_w_minus = 65;
192  u_c = 0.13;
193  u_v = 0.025;
194  u_r = 0.13;
195  u_fi = 0.13;
196  u_csi = 0.85;
197  k1 = 10;
198  k2 = 0.0;
199  break;
200  case eCF1: // (Cherry, Fenton, Am J Physiol Heart 286, 2004)
201  g_fi_max = 6.6666;
202  tau_r = 12.5;
203  tau_si = 10;
204  tau_0 = 1.5;
205  tau_v_plus = 10;
206  tau_v1_minus = 350;
207  tau_v2_minus = 80;
208  tau_w_plus = 562;
209  tau_w_minus = 48.5;
210  u_c = 0.25;
211  u_v = 0.001;
212  u_r = 0.25;
213  u_fi = 0.15;
214  u_csi = 0.2;
215  k1 = 15;
216  k2 = 0;
217  break;
218  case eCF2a: // (Cherry, Fenton, Am J Physiol Heart 286, 2004)
219  g_fi_max = 6.6666;
220  tau_r = 31;
221  tau_si = 26.5;
222  tau_0 = 1.5;
223  tau_v_plus = 10;
224  tau_v1_minus = 20;
225  tau_v2_minus = 20;
226  tau_w_plus = 800;
227  tau_w_minus = 45;
228  u_c = 0.25;
229  u_v = 0.05;
230  u_r = 0.6;
231  u_fi = 0.11;
232  u_csi = 0.7;
233  k1 = 10;
234  k2 = 1;
235  break;
236  case eCF2b: // (Cherry, Fenton, Am J Physiol Heart 286, 2004)
237  g_fi_max = 6.6666;
238  tau_r = 31;
239  tau_si = 26.5;
240  tau_0 = 1.5;
241  tau_v_plus = 10;
242  tau_v1_minus = 100;
243  tau_v2_minus = 20;
244  tau_w_plus = 800;
245  tau_w_minus = 45;
246  u_c = 0.25;
247  u_v = 0.05;
248  u_r = 0.6;
249  u_fi = 0.11;
250  u_csi = 0.7;
251  k1 = 10;
252  k2 = 1;
253  break;
254  case eCF2c: // (Cherry, Fenton, Am J Physiol Heart 286, 2004)
255  g_fi_max = 6.6666;
256  tau_r = 31;
257  tau_si = 26.5;
258  tau_0 = 1.5;
259  tau_v_plus = 10;
260  tau_v1_minus = 150;
261  tau_v2_minus = 20;
262  tau_w_plus = 800;
263  tau_w_minus = 45;
264  u_c = 0.25;
265  u_v = 0.05;
266  u_r = 0.6;
267  u_fi = 0.11;
268  u_csi = 0.7;
269  k1 = 10;
270  k2 = 1;
271  break;
272  case eCF3a: // (Cherry, Fenton, Am J Physiol Heart 286, 2004)
273  g_fi_max = 13.3333;
274  tau_r = 38;
275  tau_si = 127;
276  tau_0 = 8.3;
277  tau_v_plus = 3.33;
278  tau_v1_minus = 45;
279  tau_v2_minus = 300;
280  tau_w_plus = 600;
281  tau_w_minus = 40;
282  tau_y_plus = 1000;
283  tau_y_minus = 230;
284  u_c = 0.25;
285  u_v = 0.5;
286  u_r = 0.25;
287  u_fi = 0.25;
288  u_csi = 0.7;
289  k1 = 60;
290  k2 = 0;
291  break;
292  case eCF3b: // (Cherry, Fenton, Am J Physiol Heart 286, 2004)
293  g_fi_max = 13.3333;
294  tau_r = 38;
295  tau_si = 127;
296  tau_0 = 8.3;
297  tau_v_plus = 3.33;
298  tau_v1_minus = 20;
299  tau_v2_minus = 300;
300  tau_w_plus = 600;
301  tau_w_minus = 40;
302  tau_y_plus = 1000;
303  tau_y_minus = 230;
304  u_c = 0.25;
305  u_v = 0.5;
306  u_r = 0.25;
307  u_fi = 0.25;
308  u_csi = 0.7;
309  k1 = 60;
310  k2 = 0;
311  break;
312  case eFC2002Set1a: // (Fenton, Cherry, Chaos 12(852), 2002)
313  g_fi_max = 1/0.41; //2.44 - Fig 1A
314  tau_r = 50;
315  tau_si = 45;
316  tau_0 = 1/0.12; //8.33
317  tau_v_plus = 3.33;
318  tau_v1_minus = 1000;
319  tau_v2_minus = 19.6;
320  tau_w_plus = 667;
321  tau_w_minus = 11;
322  u_c = 0.13;
323  u_v = 0.055;
324  u_r = u_c;
325  u_fi = u_c;
326  u_csi = 0.85;
327  k1 = 10;
328  k2 = 0;
329  break;
330  case eFC2002Set1b: // (Fenton, Cherry, Chaos 12(852), 2002)
331  g_fi_max = 1/0.392; //2.55 - Fig 1B
332  tau_r = 50;
333  tau_si = 45;
334  tau_0 = 1/0.12; //8.33
335  tau_v_plus = 3.33;
336  tau_v1_minus = 1000;
337  tau_v2_minus = 19.6;
338  tau_w_plus = 667;
339  tau_w_minus = 11;
340  u_c = 0.13;
341  u_v = 0.055;
342  u_r = u_c;
343  u_fi = u_c;
344  u_csi = 0.85;
345  k1 = 10;
346  k2 = 0;
347  break;
348  case eFC2002Set1c: // (Fenton, Cherry, Chaos 12(852), 2002)
349  g_fi_max = 1/0.381; //2.62 - Fig 1C
350  tau_r = 50;
351  tau_si = 45;
352  tau_0 = 1/0.12; //8.33
353  tau_v_plus = 3.33;
354  tau_v1_minus = 1000;
355  tau_v2_minus = 19.6;
356  tau_w_plus = 667;
357  tau_w_minus = 11;
358  u_c = 0.13;
359  u_v = 0.055;
360  u_r = u_c;
361  u_fi = u_c;
362  u_csi = 0.85;
363  k1 = 10;
364  k2 = 0;
365  break;
366  case eFC2002Set1d: // (Fenton, Cherry, Chaos 12(852), 2002)
367  g_fi_max = 1/0.36; //2.77 - Fig 1D
368  tau_r = 50;
369  tau_si = 45;
370  tau_0 = 1/0.12; //8.33
371  tau_v_plus = 3.33;
372  tau_v1_minus = 1000;
373  tau_v2_minus = 19.6;
374  tau_w_plus = 667;
375  tau_w_minus = 11;
376  u_c = 0.13;
377  u_v = 0.055;
378  u_r = u_c;
379  u_fi = u_c;
380  u_csi = 0.85;
381  k1 = 10;
382  k2 = 0;
383  break;
384  case eFC2002Set1e: // (Fenton, Cherry, Chaos 12(852), 2002)
385  g_fi_max = 1/0.25; // 4 - Fig 1E
386  tau_r = 50;
387  tau_si = 45;
388  tau_0 = 1/0.12; //8.33
389  tau_v_plus = 3.33;
390  tau_v1_minus = 1000;
391  tau_v2_minus = 19.6;
392  tau_w_plus = 667;
393  tau_w_minus = 11;
394  u_c = 0.13;
395  u_v = 0.055;
396  u_r = u_c;
397  u_fi = u_c;
398  u_csi = 0.85;
399  k1 = 10;
400  k2 = 0;
401  break;
402  case eFC2002Set2: // (Fenton, Cherry, Chaos 12(852), 2002)
403  g_fi_max = 4;
404  tau_r = 190;
405  tau_si = 100000;
406  tau_0 = 10;
407  tau_v_plus = 10;
408  tau_v1_minus = 10;
409  tau_v2_minus = 10;
410  tau_w_plus = 100000;
411  tau_w_minus = 100000;
412  u_c = 0.13;
413  u_v = 100;
414  u_r = u_c;
415  u_fi = u_c;
416  u_csi = 100;
417  k1 = 10;
418  k2 = 0;
419  break;
420  case eFC2002Set4: // (Fenton, Cherry, Chaos 12(852), 2002)
421  g_fi_max = 2.46;
422  tau_r = 34;
423  tau_si = 26.5;
424  tau_0 = 9;
425  tau_v_plus = 3.33;
426  tau_v1_minus = 5; // Opposite convention
427  tau_v2_minus = 15.6; // to the paper
428  tau_w_plus = 350;
429  tau_w_minus = 80;
430  u_c = 0.15;
431  u_v = 0.04;
432  u_r = u_c;
433  u_fi = u_c;
434  u_csi = 0.45;
435  k1 = 15;
436  k2 = 0;
437  break;
438  case eFC2002Set5: // (Fenton, Cherry, Chaos 12(852), 2002)
439  g_fi_max = 2.76;
440  tau_r = 33.33;
441  tau_si = 29;
442  tau_0 = 5;
443  tau_v_plus = 3.33;
444  tau_v1_minus = 2; // Opposite convention
445  tau_v2_minus = 12; // to the paper
446  tau_w_plus = 1000;
447  tau_w_minus = 100;
448  u_c = 0.13;
449  u_v = 0.14;
450  u_r = u_c;
451  u_fi = u_c;
452  u_csi = 0.70;
453  k1 = 15;
454  k2 = 0;
455  break;
456  case eFC2002Set6: // (Fenton, Cherry, Chaos 12(852), 2002)
457  g_fi_max = 2.53;
458  tau_r = 33.33;
459  tau_si = 29;
460  tau_0 = 9;
461  tau_v_plus = 3.33;
462  tau_v1_minus = 8; // Opposite convention
463  tau_v2_minus = 9; // to the paper
464  tau_w_plus = 250;
465  tau_w_minus = 60;
466  u_c = 0.13;
467  u_v = 0.04;
468  u_r = u_c;
469  u_fi = u_c;
470  u_csi = 0.50;
471  k1 = 15;
472  k2 = 0;
473  break;
474  case eFC2002Set7: // (Fenton, Cherry, Chaos 12(852), 2002)
475  g_fi_max = 4;
476  tau_r = 100;
477  tau_si = 100000;
478  tau_0 = 12;
479  tau_v_plus = 10;
480  tau_v1_minus = 7; // Opposite convention
481  tau_v2_minus = 7; // to the paper
482  tau_w_plus = 100000;
483  tau_w_minus = 100000;
484  u_c = 0.13;
485  u_v = 100;
486  u_r = u_c;
487  u_fi = u_c;
488  u_csi = 100;
489  k1 = 10;
490  k2 = 0;
491  break;
492  case eFC2002Set8: // (Fenton, Cherry, Chaos 12(852), 2002)
493  g_fi_max = 2.2;
494  tau_r = 33.25;
495  tau_si = 29;
496  tau_0 = 12.5;
497  tau_v_plus = 13.03;
498  tau_v1_minus = 1250; // Opposite convention
499  tau_v2_minus = 19.6; // to the paper
500  tau_w_plus = 800;
501  tau_w_minus = 40;
502  u_c = 0.13;
503  u_v = 0.04;
504  u_r = u_c;
505  u_fi = u_c;
506  u_csi = 0.85;
507  k1 = 10;
508  k2 = 0;
509  break;
510  case eFC2002Set9: // (Fenton, Cherry, Chaos 12(852), 2002)
511  g_fi_max = 4;
512  tau_r = 28;
513  tau_si = 29;
514  tau_0 = 12.5;
515  tau_v_plus = 3.33;
516  tau_v1_minus = 2; // Opposite convention
517  tau_v2_minus = 15; // to the paper
518  tau_w_plus = 670;
519  tau_w_minus = 61;
520  u_c = 0.13;
521  u_v = 0.05;
522  u_r = u_c;
523  u_fi = u_c;
524  u_csi = 0.45;
525  k1 = 10;
526  k2 = 0;
527  break;
528  case eLawson: // (Lawson, Front. Physiol, 2018)
529  g_fi_max = 3;
530  tau_r = 1/0.02;
531  tau_si = 1/0.0223;
532  tau_0 = 1/0.12;
533  tau_v_plus = 3.33;
534  tau_v1_minus = 1000;
535  tau_v2_minus = 19.6;
536  tau_w_plus = 667;
537  tau_w_minus = 11;
538  u_c = 0.13;
539  u_v = 0.055;
540  u_r = u_c;
541  u_fi = u_c;
542  u_csi = 0.85;
543  k1 = 10;
544  k2 = 0;
545  break;
546  case eCAF:
547  g_fi_max = 8;
548  tau_r = 70;
549  tau_si = 114;
550  tau_0 = 32.5; //8.33
551  tau_v_plus = 5.75;
552  tau_v1_minus = 60;
553  tau_v2_minus = 82.5;
554  tau_w_plus = 300;
555  tau_w_minus = 400;
556  u_c = 0.16;
557  u_v = 0.04;
558  u_r = u_c;
559  u_fi = u_c;
560  u_csi = 0.85;
561  k1 = 10;
562  k2 = 0;
563  break;
564  }
565 
566  tau_d = C_m/g_fi_max;
567 
568  isCF3 = (model_variant == eCF3a || model_variant == eCF3b);
569 
570  // List gates and concentrations
571  m_gates.push_back(1);
572  m_gates.push_back(2);
573  m_nvar = 3;
574 
575  // Cherry-Fenton 2004 Model 3 has extra gating variable
576  if (isCF3)
577  {
578  m_gates.push_back(3);
579  m_nvar = 4;
580  }
581 
582  ASSERTL0(!isCF3, "Cherry-Fenton model 3 not implemented yet.");
583  }
584 
585 
586 
587  /**
588  *
589  */
590  FentonKarma::~FentonKarma()
591  {
592 
593  }
594 
595 
596 
597  void FentonKarma::v_Update(
598  const Array<OneD, const Array<OneD, NekDouble> >&inarray,
599  Array<OneD, Array<OneD, NekDouble> >&outarray,
600  const NekDouble time )
601  {
602  ASSERTL0(inarray.get() != outarray.get(),
603  "Must have different arrays for input and output.");
604 
605  // Variables
606  // 0 u membrane potential
607  // 1 v v gate
608  // 2 w w gate
609  int n = m_nq;
610  int i = 0;
611 
612  // Declare pointers
613  const NekDouble *u = &inarray[0][0];
614  const NekDouble *v = &inarray[1][0];
615  const NekDouble *w = &inarray[2][0];
616  const NekDouble *y = isCF3 ? &inarray[3][0] : 0;
617  NekDouble *u_new = &outarray[0][0];
618  NekDouble *v_new = &outarray[1][0];
619  NekDouble *w_new = &outarray[2][0];
620  //NekDouble *y_new = isCF3 ? &outarray[3][0] : 0;
621  NekDouble *v_tau = &m_gates_tau[0][0];
622  NekDouble *w_tau = &m_gates_tau[1][0];
623  //NekDouble *y_tau = isCF3 ? &m_gates_tau[2][0] : 0;
624 
625  // Temporary variables
626  NekDouble J_fi, J_so, J_si, h1, h2, h3, alpha, beta;
627 
628  double V;
629  // Compute rates for each point in domain
630  for (i = 0; i < n; ++i)
631  {
632  // *u is dimensional
633  // V is non-dimensional for what follows
634  V = (*u - V_0)/(V_fi - V_0);
635 
636  // Heavyside functions
637  h1 = (V < u_c) ? 0.0 : 1.0;
638  h2 = (V < u_v) ? 0.0 : 1.0;
639  h3 = (V < u_r) ? 0.0 : 1.0;
640 
641  // w-gate
642  alpha = (1-h1)/tau_w_minus;
643  beta = h1/tau_w_plus;
644  *w_tau = 1.0 / (alpha + beta);
645  *w_new = alpha * (*w_tau);
646 
647  // v-gate
648  alpha = (1-h1)/(h2*tau_v1_minus + (1-h2)*tau_v2_minus);
649  beta = h1/tau_v_plus;
650  *v_tau = 1.0 / (alpha + beta);
651  *v_new = alpha * (*v_tau);
652 
653  // y-gate
654  if (isCF3)
655  {
656  // TODO: implementation for y_tau and y_new
657  }
658 
659  // J_fi
660  J_fi = -(*v)*h1*(1 - V)*(V - u_c)/tau_d;
661 
662  // J_so
663  // added extra (1-k2*v) term from Cherry&Fenton 2004
664  J_so = V*(1-h3)*(1-k2*(*v))/tau_0 +
665  (isCF3 ? h3*V*(*y)/tau_r : h3/tau_r);
666 
667  // J_si
668  J_si = -(*w)*(1 + tanh(k1*(V - u_csi)))/(2.0*tau_si);
669 
670  // u
671  *u_new = -J_fi - J_so - J_si;
672  *u_new *= C_m*(V_fi - V_0);
673 
674  ++u, ++v, ++w, ++u_new, ++v_new, ++w_new, ++v_tau, ++w_tau;
675  }
676  }
677 
678  void FentonKarma::v_GenerateSummary(SummaryList& s)
679  {
680  SolverUtils::AddSummaryItem(s, "Cell model", "Fenton Karma");
681  SolverUtils::AddSummaryItem(s, "Cell model var.", lookupIds[model_variant]);
682  }
683 
684 
685  void FentonKarma::v_SetInitialConditions()
686  {
687  Vmath::Fill(m_nq, 0.0, m_cellSol[0], 1);
688  Vmath::Fill(m_nq, 1.0, m_cellSol[1], 1);
689  Vmath::Fill(m_nq, 1.0, m_cellSol[2], 1);
690  if (isCF3)
691  {
692  Vmath::Fill(m_nq, 0.1, m_cellSol[2], 1);
693  }
694 
695  }
696 
697  std::string FentonKarma::v_GetCellVarName(unsigned int idx)
698  {
699  switch (idx)
700  {
701  case 0: return "u";
702  case 1: return "v";
703  case 2: return "w";
704  case 3: return "y";
705  default: return "unknown";
706  }
707  }
708 }
Nektar::FentonKarma::tau_y_minus
NekDouble tau_y_minus
Definition: FentonKarma.h:96
ASSERTL0
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:216
Nektar::FentonKarma::lookupIds
static std::string lookupIds[]
Definition: FentonKarma.h:133
Nektar::CellModel::m_nq
int m_nq
Number of physical points.
Definition: CellModel.h:116
Nektar::MultiRegions::ExpListSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: LocTraceToTraceMap.h:55
Nektar::SolverUtils::SummaryList
std::vector< std::pair< std::string, std::string > > SummaryList
Definition: Misc.h:46
Vmath::Fill
void Fill(int n, const T alpha, T *x, const int incx)
Fill a vector with a constant value.
Definition: Vmath.cpp:45
Nektar::FentonKarma::~FentonKarma
virtual ~FentonKarma()
Destructor.
Definition: FentonKarma.cpp:590
VmathArray.hpp
Nektar::FentonKarma::v_GenerateSummary
virtual void v_GenerateSummary(SummaryList &s)
Prints a summary of the model parameters.
Definition: FentonKarma.cpp:678
Nektar::CellModel
Cell model base class.
Definition: CellModel.h:64
std
STL namespace.
Nektar::FentonKarma::tau_y_plus
NekDouble tau_y_plus
Definition: FentonKarma.h:95
Nektar::FentonKarma::tau_v1_minus
NekDouble tau_v1_minus
Definition: FentonKarma.h:89
Nektar::FentonKarma::v_SetInitialConditions
virtual void v_SetInitialConditions()
Definition: FentonKarma.cpp:685
Nektar::CellModel::m_gates_tau
Array< OneD, Array< OneD, NekDouble > > m_gates_tau
Storage for gate tau values.
Definition: CellModel.h:142
FentonKarma.h
Nektar::SolverUtils::AddSummaryItem
void AddSummaryItem(SummaryList &l, const std::string &name, const std::string &value)
Adds a summary item to the summary info list.
Definition: Misc.cpp:49
Nektar::CellModel::m_cellSol
Array< OneD, Array< OneD, NekDouble > > m_cellSol
Cell model solution variables.
Definition: CellModel.h:125
Nektar::CellModel::m_nvar
int m_nvar
Number of variables in cell model (inc. transmembrane voltage)
Definition: CellModel.h:118
Nektar::FentonKarma::v_Update
virtual void v_Update(const Array< OneD, const Array< OneD, NekDouble > > &inarray, Array< OneD, Array< OneD, NekDouble > > &outarray, const NekDouble time)
Computes the reaction terms $f(u,v)$ and $g(u,v)$.
Definition: FentonKarma.cpp:597
Nektar::CellModel::m_gates
std::vector< int > m_gates
Indices of cell model variables which are gates.
Definition: CellModel.h:140
Nektar::NekDouble
double NekDouble
Definition: NektarUnivTypeDefs.hpp:43
Nektar::FentonKarma::tau_w_plus
NekDouble tau_w_plus
Definition: FentonKarma.h:101
Nektar::FentonKarma::tau_v2_minus
NekDouble tau_v2_minus
Definition: FentonKarma.h:90
Nektar::GetCellModelFactory
CellModelFactory & GetCellModelFactory()
Definition: CellModel.cpp:46
Nektar::FentonKarma::model_variant
enum Variants model_variant
Definition: FentonKarma.h:131
Nektar::FentonKarma::g_fi_max
NekDouble g_fi_max
Definition: FentonKarma.h:87
Nektar::FentonKarma::tau_w_minus
NekDouble tau_w_minus
Definition: FentonKarma.h:100
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:199
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:112
Nektar::FentonKarma::v_GetCellVarName
virtual std::string v_GetCellVarName(unsigned int idx)
Definition: FentonKarma.cpp:697
Nektar::FentonKarma::tau_v_plus
NekDouble tau_v_plus
Definition: FentonKarma.h:91
Generated on Mon Dec 9 2019 21:59:37 for Nektar++ by   doxygen 1.8.13