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
40using namespace std;
41
42namespace Nektar
43{
44// Register cell model
45std::string FentonKarma::className =
46 GetCellModelFactory().RegisterCreatorFunction(
47 "FentonKarma", FentonKarma::create, "Phenomenological Model.");
48
49// Register cell model variants
50std::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
111std::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 */
122FentonKarma::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 */
649FentonKarma::~FentonKarma()
650{
651}
652
653void FentonKarma::v_Update(
654 const Array<OneD, const Array<OneD, NekDouble>> &inarray,
655 Array<OneD, Array<OneD, NekDouble>> &outarray,
656 [[maybe_unused]] const NekDouble time)
657{
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] : nullptr;
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
732void 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
738void 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
749std::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:208
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:86
Nektar::FentonKarma::tau_w_plus
NekDouble tau_w_plus
Definition: FentonKarma.h:98
Nektar::FentonKarma::tau_y_plus
NekDouble tau_y_plus
Definition: FentonKarma.h:92
Nektar::FentonKarma::v_Update
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::model_variant
enum Variants model_variant
Definition: FentonKarma.h:132
Nektar::FentonKarma::FentonKarma
FentonKarma(const LibUtilities::SessionReaderSharedPtr &pSession, const MultiRegions::ExpListSharedPtr &pField)
Constructor.
Definition: FentonKarma.cpp:122
Nektar::FentonKarma::~FentonKarma
~FentonKarma() override
Destructor.
Definition: FentonKarma.cpp:649
Nektar::FentonKarma::v_GetCellVarName
std::string v_GetCellVarName(size_t idx) override
Definition: FentonKarma.cpp:749
Nektar::FentonKarma::tau_v2_minus
NekDouble tau_v2_minus
Definition: FentonKarma.h:87
Nektar::FentonKarma::def
static std::string def
Definition: FentonKarma.h:135
Nektar::FentonKarma::create
static CellModelSharedPtr create(const LibUtilities::SessionReaderSharedPtr &pSession, const MultiRegions::ExpListSharedPtr &pField)
Creates an instance of this class.
Definition: FentonKarma.h:46
Nektar::FentonKarma::lookupIds
static std::string lookupIds[]
Definition: FentonKarma.h:134
Nektar::FentonKarma::g_fi_max
NekDouble g_fi_max
Definition: FentonKarma.h:84
Nektar::FentonKarma::tau_v_plus
NekDouble tau_v_plus
Definition: FentonKarma.h:88
Nektar::FentonKarma::v_GenerateSummary
void v_GenerateSummary(SummaryList &s) override
Prints a summary of the model parameters.
Definition: FentonKarma.cpp:732
Nektar::FentonKarma::v_SetInitialConditions
void v_SetInitialConditions() override
Definition: FentonKarma.cpp:738
Nektar::FentonKarma::tau_w_minus
NekDouble tau_w_minus
Definition: FentonKarma.h:97
Nektar::FentonKarma::tau_y_minus
NekDouble tau_y_minus
Definition: FentonKarma.h:93
Nektar::FentonKarma::className
static std::string className
Name of class.
Definition: FentonKarma.h:54
Nektar::LibUtilities::NekFactory::RegisterCreatorFunction
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, std::string pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:197
Nektar::LibUtilities::SessionReader::RegisterEnumValue
static std::string RegisterEnumValue(std::string pEnum, std::string pString, int pEnumValue)
Registers an enumeration value.
Definition: SessionReader.h:559
Nektar::LibUtilities::SessionReader::RegisterDefaultSolverInfo
static std::string RegisterDefaultSolverInfo(const std::string &pName, const std::string &pValue)
Registers the default string value of a solver info property.
Definition: SessionReader.h:593
Nektar::LibUtilities::SessionReaderSharedPtr
std::shared_ptr< SessionReader > SessionReaderSharedPtr
Definition: SessionReader.h:113
Nektar::LibUtilities::beta
@ beta
Gauss Radau pinned at x=-1,.
Definition: PointsType.h:59
Nektar::MultiRegions::ExpListSharedPtr
std::shared_ptr< ExpList > ExpListSharedPtr
Shared pointer to an ExpList object.
Definition: LocTraceToTraceMap.h:56
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::UnitTests::w
std::vector< double > w(NPUPPER)
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.hpp:54
Generated on Wed Jan 31 2024 09:35:57 for Nektar++ by doxygen 1.9.4