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