Skip to content

Commit 28be387

Browse files
javijv4ktbolt
javijv4
and
ktbolt
authored
Fixing domain-specific conductances for EP (#506)
* adding domain specific conductances * adding test for domain * adding test to test_cep.py * removing unused function parameter; clarifying comments; fixing bug in cep_ion.cpp * no domains in equations fixed; updated parameter assignment --------- Co-authored-by: Dave Parker <dvdpax@gmail.com>
1 parent bc36c39 commit 28be387

14 files changed

Lines changed: 3576 additions & 56 deletions

File tree

‎Code/Source/solver/CepMod.h‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -165,6 +165,18 @@ class cepModelType
165165

166166
/// @brief Time integration options
167167
odeType odes;
168+
169+
/// @brief Interface for Aliev-Panfilov cellular activation model
170+
CepModAp ap;
171+
172+
/// @brief Interface for Bueno-Orovio cellular activation model
173+
CepModBo bo;
174+
175+
/// @brief Interface for Fitzhugh-Nagumo cellular activation model
176+
CepModFn fn;
177+
178+
/// @brief Interface for Tusscher-Panfilov cellular activation model
179+
CepModTtp ttp;
168180
};
169181

170182
/// @brief Cardiac electromechanics model type

‎Code/Source/solver/cep_ion.cpp‎

Lines changed: 37 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -52,7 +52,7 @@ void cep_init(Simulation* simulation)
5252
for (int iDmn = 0; iDmn < eq.nDmn; iDmn++) {
5353
auto cPhys = eq.dmn[iDmn].phys;
5454
int dID = eq.dmn[iDmn].Id;
55-
if ((cPhys != EquationType::phys_CEP) || !utils::btest(com_mod.dmnId(a),dID)) {
55+
if ((cPhys != EquationType::phys_CEP) || (dID >= 0 && !utils::btest(com_mod.dmnId(a),dID))) {
5656
continue;
5757
}
5858
int nX = eq.dmn[iDmn].cep.nX;
@@ -62,7 +62,7 @@ void cep_init(Simulation* simulation)
6262
Vector<double> Xl(nX);
6363
Vector<double> Xgl(nG);
6464

65-
cep_init_l(cep_mod, eq.dmn[iDmn].cep, nX, nG, Xl, Xgl);
65+
cep_init_l(eq.dmn[iDmn].cep, nX, nG, Xl, Xgl);
6666

6767
sA(a) = sA(a) + 1.0;
6868

@@ -97,7 +97,7 @@ void cep_init(Simulation* simulation)
9797
Vector<double> Xl(nX);
9898
Vector<double> Xgl(nG);
9999

100-
cep_init_l(cep_mod, eq.dmn[1].cep, nX, nG, Xl, Xgl);
100+
cep_init_l(eq.dmn[0].cep, nX, nG, Xl, Xgl);
101101

102102
for (int i = 0; i < nX; i++) {
103103
cep_mod.Xion(i,a) = Xl(i);
@@ -114,24 +114,24 @@ void cep_init(Simulation* simulation)
114114
// cep_init_l
115115
//------------
116116
//
117-
void cep_init_l(CepMod& cep_mod, cepModelType& cep, int nX, int nG, Vector<double>& X, Vector<double>& Xg)
117+
void cep_init_l(cepModelType& cep, int nX, int nG, Vector<double>& X, Vector<double>& Xg)
118118
{
119119
switch (cep.cepType) {
120120

121121
case ElectrophysiologyModelType::AP:
122-
cep_mod.ap.init(nX, X);
122+
cep.ap.init(nX, X);
123123
break;
124124

125125
case ElectrophysiologyModelType::BO:
126-
cep_mod.bo.init(nX, X);
126+
cep.bo.init(nX, X);
127127
break;
128128

129129
case ElectrophysiologyModelType::FN:
130-
cep_mod.fn.init(nX, X);
130+
cep.fn.init(nX, X);
131131
break;
132132

133133
case ElectrophysiologyModelType::TTP:
134-
cep_mod.ttp.init(cep.imyo, nX, nG, X, Xg);
134+
cep.ttp.init(cep.imyo, nX, nG, X, Xg);
135135
break;
136136
}
137137
}
@@ -220,9 +220,9 @@ void cep_integ(Simulation* simulation, const int iEq, const int iDof, const Arra
220220
auto cPhys = dmn.phys;
221221
int dID = dmn.Id;
222222

223-
if (cPhys != Equation_CEP || !utils::btest(com_mod.dmnId(Ac),dID)) {
223+
if (cPhys != Equation_CEP || (dID >= 0 && !utils::btest(com_mod.dmnId(Ac),dID))) {
224224
continue;
225-
}
225+
}
226226

227227
int nX = dmn.cep.nX;
228228
int nG = dmn.cep.nG;
@@ -383,12 +383,12 @@ void cep_integ_l(CepMod& cep_mod, cepModelType& cep, int nX, int nG, Vector<doub
383383
} else {
384384
Istim = 0.0;
385385
}
386-
cep_mod.ap.integ_fe(nX, X, t, cep.dt, Istim, Ksac);
386+
cep.ap.integ_fe(nX, X, t, cep.dt, Istim, Ksac);
387387

388388
// Electromechanics excitation-activation
389389
if (cem.aStress) {
390390
double epsX;
391-
cep_mod.ap.actv_strs(X(0), cep.dt, yl, epsX);
391+
cep.ap.actv_strs(X(0), cep.dt, yl, epsX);
392392
}
393393
}
394394
} break;
@@ -402,12 +402,12 @@ void cep_integ_l(CepMod& cep_mod, cepModelType& cep, int nX, int nG, Vector<doub
402402
} else {
403403
Istim = 0.0;
404404
}
405-
cep_mod.ap.integ_rk(nX, X, t, cep.dt, Istim, Ksac);
405+
cep.ap.integ_rk(nX, X, t, cep.dt, Istim, Ksac);
406406

407407
// Electromechanics excitation-activation
408408
if (cem.aStress) {
409409
double epsX;
410-
cep_mod.ap.actv_strs(X(0), cep.dt, yl, epsX);
410+
cep.ap.actv_strs(X(0), cep.dt, yl, epsX);
411411
}
412412
}
413413
} break;
@@ -422,12 +422,12 @@ void cep_integ_l(CepMod& cep_mod, cepModelType& cep, int nX, int nG, Vector<doub
422422
Istim = 0.0;
423423
}
424424

425-
cep_mod.ap.integ_cn2(nX, X, t, cep.dt, Istim, Ksac, IPAR, RPAR);
425+
cep.ap.integ_cn2(nX, X, t, cep.dt, Istim, Ksac, IPAR, RPAR);
426426

427427
// Electromechanics excitation-activation
428428
if (cem.aStress) {
429429
double epsX;
430-
cep_mod.ap.actv_strs(X(0), cep.dt, yl, epsX);
430+
cep.ap.actv_strs(X(0), cep.dt, yl, epsX);
431431
}
432432
}
433433
} break;
@@ -453,14 +453,14 @@ void cep_integ_l(CepMod& cep_mod, cepModelType& cep, int nX, int nG, Vector<doub
453453
Istim = 0.0;
454454
}
455455

456-
cep_mod.bo.integ_fe(cep.imyo, nX, X, t, cep.dt, Istim, Ksac, RPAR);
456+
cep.bo.integ_fe(cep.imyo, nX, X, t, cep.dt, Istim, Ksac, RPAR);
457457

458458
// Electromechanics excitation-activation
459459
if (cem.aStress) {
460460
double epsX;
461-
cep_mod.bo.actv_strs(X(0), cep.dt, yl, epsX);
461+
cep.bo.actv_strs(X(0), cep.dt, yl, epsX);
462462
} else if (cem.aStrain) {
463-
cep_mod.bo.actv_strn(X(3), I4f, cep.dt, yl);
463+
cep.bo.actv_strn(X(3), I4f, cep.dt, yl);
464464
}
465465
}
466466
} break;
@@ -475,14 +475,14 @@ void cep_integ_l(CepMod& cep_mod, cepModelType& cep, int nX, int nG, Vector<doub
475475
Istim = 0.0;
476476
}
477477

478-
cep_mod.bo.integ_rk(cep.imyo, nX, X, t, cep.dt, Istim, Ksac, RPAR);
478+
cep.bo.integ_rk(cep.imyo, nX, X, t, cep.dt, Istim, Ksac, RPAR);
479479

480480
// Electromechanics excitation-activation
481481
if (cem.aStress) {
482482
double epsX;
483-
cep_mod.bo.actv_strs(X(0), cep.dt, yl, epsX);
483+
cep.bo.actv_strs(X(0), cep.dt, yl, epsX);
484484
} else if (cem.aStrain) {
485-
cep_mod.bo.actv_strn(X(3), I4f, cep.dt, yl);
485+
cep.bo.actv_strn(X(3), I4f, cep.dt, yl);
486486
}
487487
}
488488
} break;
@@ -497,14 +497,14 @@ void cep_integ_l(CepMod& cep_mod, cepModelType& cep, int nX, int nG, Vector<doub
497497
Istim = 0.0;
498498
}
499499

500-
cep_mod.bo.integ_cn2(cep.imyo, nX, X, t, cep.dt, Istim, Ksac, IPAR, RPAR);
500+
cep.bo.integ_cn2(cep.imyo, nX, X, t, cep.dt, Istim, Ksac, IPAR, RPAR);
501501

502502
// Electromechanics excitation-activation
503503
if (cem.aStress) {
504504
double epsX;
505-
cep_mod.bo.actv_strs(X(0), cep.dt, yl, epsX);
505+
cep.bo.actv_strs(X(0), cep.dt, yl, epsX);
506506
} else if (cem.aStrain) {
507-
cep_mod.bo.actv_strn(X(3), I4f, cep.dt, yl);
507+
cep.bo.actv_strn(X(3), I4f, cep.dt, yl);
508508
}
509509
}
510510
} break;
@@ -530,7 +530,7 @@ void cep_integ_l(CepMod& cep_mod, cepModelType& cep, int nX, int nG, Vector<doub
530530
} else {
531531
Istim = 0.0;
532532
}
533-
cep_mod.fn.integ_fe(nX, X, t, cep.dt, Istim);
533+
cep.fn.integ_fe(nX, X, t, cep.dt, Istim);
534534
}
535535
} break;
536536

@@ -543,7 +543,7 @@ void cep_integ_l(CepMod& cep_mod, cepModelType& cep, int nX, int nG, Vector<doub
543543
} else {
544544
Istim = 0.0;
545545
}
546-
cep_mod.fn.integ_rk(nX, X, t, cep.dt, Istim);
546+
cep.fn.integ_rk(nX, X, t, cep.dt, Istim);
547547
}
548548
} break;
549549

@@ -556,7 +556,7 @@ void cep_integ_l(CepMod& cep_mod, cepModelType& cep, int nX, int nG, Vector<doub
556556
} else {
557557
Istim = 0.0;
558558
}
559-
cep_mod.fn.integ_cn2(nX, X, t, cep.dt, Istim, IPAR, RPAR);
559+
cep.fn.integ_cn2(nX, X, t, cep.dt, Istim, IPAR, RPAR);
560560
}
561561
} break;
562562
}
@@ -580,14 +580,14 @@ void cep_integ_l(CepMod& cep_mod, cepModelType& cep, int nX, int nG, Vector<doub
580580
} else {
581581
Istim = 0.0;
582582
}
583-
cep_mod.ttp.integ_fe(cep.imyo, nX, nG, X, Xg, t, cep.dt, Istim, Ksac, RPAR);
583+
cep.ttp.integ_fe(cep.imyo, nX, nG, X, Xg, t, cep.dt, Istim, Ksac, RPAR);
584584

585585
// Electromechanics excitation-activation
586586
if (cem.aStress) {
587587
double epsX;
588-
cep_mod.ttp.actv_strs(X(3), cep.dt, yl, epsX);
588+
cep.ttp.actv_strs(X(3), cep.dt, yl, epsX);
589589
} else if (cem.aStrain) {
590-
cep_mod.ttp.actv_strn(X(3), I4f, cep.dt, yl);
590+
cep.ttp.actv_strn(X(3), I4f, cep.dt, yl);
591591
}
592592
}
593593
} break;
@@ -602,14 +602,14 @@ void cep_integ_l(CepMod& cep_mod, cepModelType& cep, int nX, int nG, Vector<doub
602602
Istim = 0.0;
603603
}
604604

605-
cep_mod.ttp.integ_rk(cep.imyo, nX, nG, X, Xg, t, cep.dt, Istim, Ksac, RPAR);
605+
cep.ttp.integ_rk(cep.imyo, nX, nG, X, Xg, t, cep.dt, Istim, Ksac, RPAR);
606606

607607
// Electromechanics excitation-activation
608608
if (cem.aStress) {
609609
double epsX;
610-
cep_mod.ttp.actv_strs(X(3), cep.dt, yl, epsX);
610+
cep.ttp.actv_strs(X(3), cep.dt, yl, epsX);
611611
} else if (cem.aStrain) {
612-
cep_mod.ttp.actv_strn(X(3), I4f, cep.dt, yl);
612+
cep.ttp.actv_strn(X(3), I4f, cep.dt, yl);
613613
}
614614
}
615615
} break;
@@ -624,14 +624,14 @@ void cep_integ_l(CepMod& cep_mod, cepModelType& cep, int nX, int nG, Vector<doub
624624
Istim = 0.0;
625625
}
626626

627-
cep_mod.ttp.integ_cn2(cep.imyo, nX, nG, X, Xg, t, cep.dt, Istim, Ksac, IPAR, RPAR);
627+
cep.ttp.integ_cn2(cep.imyo, nX, nG, X, Xg, t, cep.dt, Istim, Ksac, IPAR, RPAR);
628628

629629
// Electromechanics excitation-activation
630630
if (cem.aStress) {
631631
double epsX;
632-
cep_mod.ttp.actv_strs(X(3), cep.dt, yl, epsX);
632+
cep.ttp.actv_strs(X(3), cep.dt, yl, epsX);
633633
} else if (cem.aStrain) {
634-
cep_mod.ttp.actv_strn(X(3), I4f, cep.dt, yl);
634+
cep.ttp.actv_strn(X(3), I4f, cep.dt, yl);
635635
}
636636
}
637637
} break;

‎Code/Source/solver/cep_ion.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -17,7 +17,7 @@ namespace cep_ion {
1717

1818
void cep_init(Simulation* simulation);
1919

20-
void cep_init_l(CepMod& cep_mod, cepModelType& cep, int nX, int nG, Vector<double>& X, Vector<double>& Xg);
20+
void cep_init_l(cepModelType& cep, int nX, int nG, Vector<double>& X, Vector<double>& Xg);
2121

2222
void cep_integ(Simulation* simulation, const int iEq, const int iDof, const Array<double>& Dg);
2323

‎Code/Source/solver/distribute.cpp‎

Lines changed: 9 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1558,14 +1558,15 @@ void dist_eq(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, const std::
15581558
cm.bcast(cm_mod, &cep.odes.relTol);
15591559
}
15601560

1561-
cm.bcast(cm_mod, &cep_mod.ttp.G_Na);
1562-
cm.bcast(cm_mod, &cep_mod.ttp.G_CaL);
1563-
cm.bcast(cm_mod, &cep_mod.ttp.G_Kr);
1564-
cm.bcast(cm_mod, cep_mod.ttp.G_Ks);
1565-
cm.bcast(cm_mod, cep_mod.ttp.G_to);
1566-
1567-
cm.bcast(cm_mod, cep_mod.bo.tau_si);
1568-
cm.bcast(cm_mod, cep_mod.bo.tau_fi);
1561+
// Broadcast domain-specific model parameters
1562+
cm.bcast(cm_mod, &cep.ttp.G_Na);
1563+
cm.bcast(cm_mod, &cep.ttp.G_CaL);
1564+
cm.bcast(cm_mod, &cep.ttp.G_Kr);
1565+
cm.bcast(cm_mod, cep.ttp.G_Ks);
1566+
cm.bcast(cm_mod, cep.ttp.G_to);
1567+
1568+
cm.bcast(cm_mod, cep.bo.tau_si);
1569+
cm.bcast(cm_mod, cep.bo.tau_fi);
15691570
}
15701571

15711572
if ((dmn.phys == EquationType::phys_struct) || (dmn.phys == EquationType::phys_ustruct)) {

‎Code/Source/solver/read_files.cpp‎

Lines changed: 19 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -1000,17 +1000,26 @@ void read_cep_domain(Simulation* simulation, EquationParameters* eq_params, Doma
10001000
}
10011001

10021002
// Set Ttp parameters.
1003-
//
1004-
if (domain_params->G_Na.defined()) { cep_mod.ttp.G_Na = domain_params->G_Na.value(); }
1005-
if (domain_params->G_Kr.defined()) { cep_mod.ttp.G_Kr = domain_params->G_Kr.value(); }
1006-
if (domain_params->G_Ks.defined()) { cep_mod.ttp.G_Ks[lDmn.cep.imyo - 1] = domain_params->G_Ks.value(); }
1007-
if (domain_params->G_to.defined()) { cep_mod.ttp.G_to[lDmn.cep.imyo - 1] = domain_params->G_to.value(); }
1008-
if (domain_params->G_CaL.defined()) { cep_mod.ttp.G_CaL = domain_params->G_CaL.value(); }
1003+
// Scalar params
1004+
std::map<Parameter<double>*,double*> simple_ttp_params{
1005+
{&domain_params->G_Na, &lDmn.cep.ttp.G_Na},
1006+
{&domain_params->G_Kr, &lDmn.cep.ttp.G_Kr},
1007+
{&domain_params->G_CaL, &lDmn.cep.ttp.G_CaL}
1008+
};
10091009

1010-
// Set Bo parameters.
1011-
//
1012-
if (domain_params->tau_si.defined()) { cep_mod.bo.tau_si[lDmn.cep.imyo - 1] = domain_params->tau_si.value(); }
1013-
if (domain_params->tau_fi.defined()) { cep_mod.bo.tau_fi[lDmn.cep.imyo - 1] = domain_params->tau_fi.value(); }
1010+
for (auto& [param, value] : simple_ttp_params) {
1011+
if (param->defined()) {
1012+
*value = param->value();
1013+
}
1014+
}
1015+
1016+
// Array params
1017+
if (domain_params->G_Ks.defined()) {
1018+
lDmn.cep.ttp.G_Ks[lDmn.cep.imyo - 1] = domain_params->G_Ks.value();
1019+
}
1020+
if (domain_params->G_to.defined()) {
1021+
lDmn.cep.ttp.G_to[lDmn.cep.imyo - 1] = domain_params->G_to.value();
1022+
}
10141023

10151024
// Set stimulus parameters.
10161025
//

‎tests/cases/cep/slab_domains/README.md‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,11 @@
1+
2+
# **Problem Description**
3+
4+
Simulate the propagation of a planar wave in a slab with different domain conductances using the ten-Tusscher-Panfilov cell activation model.
5+
6+
![Slab domains visualization](domains.png)
7+
8+
## References
9+
K. H. W. J. ten Tusscher, D. Noble, P. J. Noble, and A. V. Panfilov. A model for human ventricular tissue. American Journal of Physiology-Heart and Circulatory Physiology,
10+
286(4):H1573–H1589, apr 2004.
11+

‎tests/cases/cep/slab_domains/domains.png‎

Lines changed: 3 additions & 0 deletions
Loading

‎tests/cases/cep/slab_domains/mesh/VOLUME.vtu‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,3 @@
1+
version https://git-lfs.github.com/spec/v1
2+
oid sha256:69bb6a21296fa79919d0a19669a33a28ae28efdb1bd38124f854815d5d406e8f
3+
size 89215

0 commit comments

Comments
 (0)