Skip to content

Fixing domain-specific conductances for EP #506

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
javijv4 wants to merge 6 commits into
base: main
Choose a base branch
from
Open

Fixing domain-specific conductances for EP #506

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
00f7853
adding domain specific conductances
javijv4 Feb 11, 2026
f46eedc
adding test for domain
javijv4 Feb 11, 2026
cf07869
Merge branch 'SimVascular:main' into dev_#504
javijv4 Feb 11, 2026
15c866d
adding test to test_cep.py
javijv4 Feb 12, 2026
06a3f33
removing unused function parameter; clarifying comments; fixing bug i…
javijv4 Feb 12, 2026
bb4887b
no domains in equations fixed; updated parameter assignment
javijv4 Feb 19, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 12 additions & 0 deletions Code/Source/solver/CepMod.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -165,6 +165,18 @@ class cepModelType

/// @brief Time integration options
odeType odes;

/// @brief Domain-specific interface for Aliev-Panfilov cellular activation model
CepModAp ap;

/// @brief Domain-specific interface for Bueno-Orovio cellular activation model
CepModBo bo;

/// @brief Domain-specific interface for Fitzhugh-Nagumo cellular activation model
CepModFn fn;

/// @brief Domain-specific interface for Tusscher-Panfilov cellular activation model
CepModTtp ttp;
};

/// @brief Cardiac electromechanics model type
Expand Down
62 changes: 31 additions & 31 deletions Code/Source/solver/cep_ion.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -119,19 +119,19 @@ void cep_init_l(CepMod& cep_mod, cepModelType& cep, int nX, int nG, Vector<doubl
switch (cep.cepType) {

case ElectrophysiologyModelType::AP:
cep_mod.ap.init(nX, X);
cep.ap.init(nX, X);
break;

case ElectrophysiologyModelType::BO:
cep_mod.bo.init(nX, X);
cep.bo.init(nX, X);
break;

case ElectrophysiologyModelType::FN:
cep_mod.fn.init(nX, X);
cep.fn.init(nX, X);
break;

case ElectrophysiologyModelType::TTP:
cep_mod.ttp.init(cep.imyo, nX, nG, X, Xg);
cep.ttp.init(cep.imyo, nX, nG, X, Xg);
break;
}
}
Expand Down Expand Up @@ -383,12 +383,12 @@ void cep_integ_l(CepMod& cep_mod, cepModelType& cep, int nX, int nG, Vector<doub
} else {
Istim = 0.0;
}
cep_mod.ap.integ_fe(nX, X, t, cep.dt, Istim, Ksac);
cep.ap.integ_fe(nX, X, t, cep.dt, Istim, Ksac);

// Electromechanics excitation-activation
if (cem.aStress) {
double epsX;
cep_mod.ap.actv_strs(X(0), cep.dt, yl, epsX);
cep.ap.actv_strs(X(0), cep.dt, yl, epsX);
}
}
} break;
Expand All @@ -402,12 +402,12 @@ void cep_integ_l(CepMod& cep_mod, cepModelType& cep, int nX, int nG, Vector<doub
} else {
Istim = 0.0;
}
cep_mod.ap.integ_rk(nX, X, t, cep.dt, Istim, Ksac);
cep.ap.integ_rk(nX, X, t, cep.dt, Istim, Ksac);

// Electromechanics excitation-activation
if (cem.aStress) {
double epsX;
cep_mod.ap.actv_strs(X(0), cep.dt, yl, epsX);
cep.ap.actv_strs(X(0), cep.dt, yl, epsX);
}
}
} break;
Expand All @@ -422,12 +422,12 @@ void cep_integ_l(CepMod& cep_mod, cepModelType& cep, int nX, int nG, Vector<doub
Istim = 0.0;
}

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

// Electromechanics excitation-activation
if (cem.aStress) {
double epsX;
cep_mod.ap.actv_strs(X(0), cep.dt, yl, epsX);
cep.ap.actv_strs(X(0), cep.dt, yl, epsX);
}
}
} break;
Expand All @@ -453,14 +453,14 @@ void cep_integ_l(CepMod& cep_mod, cepModelType& cep, int nX, int nG, Vector<doub
Istim = 0.0;
}

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

// Electromechanics excitation-activation
if (cem.aStress) {
double epsX;
cep_mod.bo.actv_strs(X(0), cep.dt, yl, epsX);
cep.bo.actv_strs(X(0), cep.dt, yl, epsX);
} else if (cem.aStrain) {
cep_mod.bo.actv_strn(X(3), I4f, cep.dt, yl);
cep.bo.actv_strn(X(3), I4f, cep.dt, yl);
}
}
} break;
Expand All @@ -475,14 +475,14 @@ void cep_integ_l(CepMod& cep_mod, cepModelType& cep, int nX, int nG, Vector<doub
Istim = 0.0;
}

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

// Electromechanics excitation-activation
if (cem.aStress) {
double epsX;
cep_mod.bo.actv_strs(X(0), cep.dt, yl, epsX);
cep.bo.actv_strs(X(0), cep.dt, yl, epsX);
} else if (cem.aStrain) {
cep_mod.bo.actv_strn(X(3), I4f, cep.dt, yl);
cep.bo.actv_strn(X(3), I4f, cep.dt, yl);
}
}
} break;
Expand All @@ -497,14 +497,14 @@ void cep_integ_l(CepMod& cep_mod, cepModelType& cep, int nX, int nG, Vector<doub
Istim = 0.0;
}

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

// Electromechanics excitation-activation
if (cem.aStress) {
double epsX;
cep_mod.bo.actv_strs(X(0), cep.dt, yl, epsX);
cep.bo.actv_strs(X(0), cep.dt, yl, epsX);
} else if (cem.aStrain) {
cep_mod.bo.actv_strn(X(3), I4f, cep.dt, yl);
cep.bo.actv_strn(X(3), I4f, cep.dt, yl);
}
}
} break;
Expand All @@ -530,7 +530,7 @@ void cep_integ_l(CepMod& cep_mod, cepModelType& cep, int nX, int nG, Vector<doub
} else {
Istim = 0.0;
}
cep_mod.fn.integ_fe(nX, X, t, cep.dt, Istim);
cep.fn.integ_fe(nX, X, t, cep.dt, Istim);
}
} break;

Expand All @@ -543,7 +543,7 @@ void cep_integ_l(CepMod& cep_mod, cepModelType& cep, int nX, int nG, Vector<doub
} else {
Istim = 0.0;
}
cep_mod.fn.integ_rk(nX, X, t, cep.dt, Istim);
cep.fn.integ_rk(nX, X, t, cep.dt, Istim);
}
} break;

Expand All @@ -556,7 +556,7 @@ void cep_integ_l(CepMod& cep_mod, cepModelType& cep, int nX, int nG, Vector<doub
} else {
Istim = 0.0;
}
cep_mod.fn.integ_cn2(nX, X, t, cep.dt, Istim, IPAR, RPAR);
cep.fn.integ_cn2(nX, X, t, cep.dt, Istim, IPAR, RPAR);
}
} break;
}
Expand All @@ -580,14 +580,14 @@ void cep_integ_l(CepMod& cep_mod, cepModelType& cep, int nX, int nG, Vector<doub
} else {
Istim = 0.0;
}
cep_mod.ttp.integ_fe(cep.imyo, nX, nG, X, Xg, t, cep.dt, Istim, Ksac, RPAR);
cep.ttp.integ_fe(cep.imyo, nX, nG, X, Xg, t, cep.dt, Istim, Ksac, RPAR);

// Electromechanics excitation-activation
if (cem.aStress) {
double epsX;
cep_mod.ttp.actv_strs(X(3), cep.dt, yl, epsX);
cep.ttp.actv_strs(X(3), cep.dt, yl, epsX);
} else if (cem.aStrain) {
cep_mod.ttp.actv_strn(X(3), I4f, cep.dt, yl);
cep.ttp.actv_strn(X(3), I4f, cep.dt, yl);
}
}
} break;
Expand All @@ -602,14 +602,14 @@ void cep_integ_l(CepMod& cep_mod, cepModelType& cep, int nX, int nG, Vector<doub
Istim = 0.0;
}

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

// Electromechanics excitation-activation
if (cem.aStress) {
double epsX;
cep_mod.ttp.actv_strs(X(3), cep.dt, yl, epsX);
cep.ttp.actv_strs(X(3), cep.dt, yl, epsX);
} else if (cem.aStrain) {
cep_mod.ttp.actv_strn(X(3), I4f, cep.dt, yl);
cep.ttp.actv_strn(X(3), I4f, cep.dt, yl);
}
}
} break;
Expand All @@ -624,14 +624,14 @@ void cep_integ_l(CepMod& cep_mod, cepModelType& cep, int nX, int nG, Vector<doub
Istim = 0.0;
}

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

// Electromechanics excitation-activation
if (cem.aStress) {
double epsX;
cep_mod.ttp.actv_strs(X(3), cep.dt, yl, epsX);
cep.ttp.actv_strs(X(3), cep.dt, yl, epsX);
} else if (cem.aStrain) {
cep_mod.ttp.actv_strn(X(3), I4f, cep.dt, yl);
cep.ttp.actv_strn(X(3), I4f, cep.dt, yl);
}
}
} break;
Expand Down
17 changes: 9 additions & 8 deletions Code/Source/solver/distribute.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1557,14 +1557,15 @@ void dist_eq(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, const std::
cm.bcast(cm_mod, &cep.odes.relTol);
}

cm.bcast(cm_mod, &cep_mod.ttp.G_Na);
cm.bcast(cm_mod, &cep_mod.ttp.G_CaL);
cm.bcast(cm_mod, &cep_mod.ttp.G_Kr);
cm.bcast(cm_mod, cep_mod.ttp.G_Ks);
cm.bcast(cm_mod, cep_mod.ttp.G_to);

cm.bcast(cm_mod, cep_mod.bo.tau_si);
cm.bcast(cm_mod, cep_mod.bo.tau_fi);
// Broadcast domain-specific model parameters
cm.bcast(cm_mod, &cep.ttp.G_Na);
cm.bcast(cm_mod, &cep.ttp.G_CaL);
cm.bcast(cm_mod, &cep.ttp.G_Kr);
cm.bcast(cm_mod, cep.ttp.G_Ks);
cm.bcast(cm_mod, cep.ttp.G_to);

cm.bcast(cm_mod, cep.bo.tau_si);
cm.bcast(cm_mod, cep.bo.tau_fi);
}

if ((dmn.phys == EquationType::phys_struct) || (dmn.phys == EquationType::phys_ustruct)) {
Expand Down
18 changes: 9 additions & 9 deletions Code/Source/solver/read_files.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -999,18 +999,18 @@ void read_cep_domain(Simulation* simulation, EquationParameters* eq_params, Doma
}
}

// Set Ttp parameters.
// Set Ttp parameters (domain-specific).
//
if (domain_params->G_Na.defined()) { cep_mod.ttp.G_Na = domain_params->G_Na.value(); }
if (domain_params->G_Kr.defined()) { cep_mod.ttp.G_Kr = domain_params->G_Kr.value(); }
if (domain_params->G_Ks.defined()) { cep_mod.ttp.G_Ks[lDmn.cep.imyo - 1] = domain_params->G_Ks.value(); }
if (domain_params->G_to.defined()) { cep_mod.ttp.G_to[lDmn.cep.imyo - 1] = domain_params->G_to.value(); }
if (domain_params->G_CaL.defined()) { cep_mod.ttp.G_CaL = domain_params->G_CaL.value(); }
if (domain_params->G_Na.defined()) { lDmn.cep.ttp.G_Na = domain_params->G_Na.value(); }
if (domain_params->G_Kr.defined()) { lDmn.cep.ttp.G_Kr = domain_params->G_Kr.value(); }
if (domain_params->G_Ks.defined()) { lDmn.cep.ttp.G_Ks[lDmn.cep.imyo - 1] = domain_params->G_Ks.value(); }
if (domain_params->G_to.defined()) { lDmn.cep.ttp.G_to[lDmn.cep.imyo - 1] = domain_params->G_to.value(); }
Copy link
Collaborator

@ktbolt ktbolt Feb 13, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@javijv4 I don't like the look of this

lDmn.cep.ttp.G_to[lDmn.cep.imyo - 1]

cepModelType::imyo should not even be an int, should be an enum.

Copy link
Contributor Author

@javijv4 javijv4 Feb 17, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ktbolt I also don't like it. I believe this statement lDmn.cep.ttp.G_to[lDmn.cep.imyo - 1] is there to handle within the code different myocardial zones. I think the user should define different domains and assign parameters to each domain in the .xml file, such that there is a single G_Ks, G_to for each domain, same as with G_Na, G_Kr and G_Cal.

Can we create a separate issue for this? I think CepModTtp shouldn't use imyo at all, and the user should create separate domains to assign parameters to different regions, with some useful default values defined in the code (as in this issue here). But this is more of a TTP issue, different from the domain-specific values.

@kko27, does modifying TTP to be agnostic of imyo make sense to you? (If we can correctly assign conductances based on domains).

Copy link
Contributor

@kko27 kko27 Feb 18, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@javijv4 and @ktbolt. This was my original plan in my issue but I was thinking for backward compatibility, it made sense to keep these hard-coded values. But it works for me to remove the initial conditions entirely from the codebase

Copy link
Contributor Author

@javijv4 javijv4 Feb 19, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@kko27, we'll leave it as it is for now, and probably address it when we refactor this code.

if (domain_params->G_CaL.defined()) { lDmn.cep.ttp.G_CaL = domain_params->G_CaL.value(); }

// Set Bo parameters.
// Set Bo parameters (domain-specific).
//
if (domain_params->tau_si.defined()) { cep_mod.bo.tau_si[lDmn.cep.imyo - 1] = domain_params->tau_si.value(); }
if (domain_params->tau_fi.defined()) { cep_mod.bo.tau_fi[lDmn.cep.imyo - 1] = domain_params->tau_fi.value(); }
if (domain_params->tau_si.defined()) { lDmn.cep.bo.tau_si[lDmn.cep.imyo - 1] = domain_params->tau_si.value(); }
if (domain_params->tau_fi.defined()) { lDmn.cep.bo.tau_fi[lDmn.cep.imyo - 1] = domain_params->tau_fi.value(); }

// Set stimulus parameters.
//
Expand Down
11 changes: 11 additions & 0 deletions tests/cases/cep/slab_domains/README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,11 @@

# **Problem Description**

Simulate the propagation of a planar wave in a slab with different domain conductances using the ten-Tusscher-Panfilov cell activation model.

![Slab domains visualization](domains.png)

## References
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,
286(4):H1573–H1589, apr 2004.

3 changes: 3 additions & 0 deletions tests/cases/cep/slab_domains/domains.png
View file
Open in desktop
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
3 changes: 3 additions & 0 deletions tests/cases/cep/slab_domains/mesh/VOLUME.vtu
View file
Open in desktop
Git LFS file not shown
Loading
Loading