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