Skip to content

Deleting unused code in TTP ionic model #568

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

Merged
javijv4 merged 1 commit into from
Jun 18, 2026
Merged

Deleting unused code in TTP ionic model #568

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
241 changes: 0 additions & 241 deletions Code/Source/solver/ionic_ttp.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -456,246 +456,5 @@ Vector<double> TTP::getf(const unsigned int zone_id, const Vector<double> &X,
return dX;
}

Array<double> TTP::getj(const unsigned int zone_id, const Vector<double> &X,
const Vector<double> &Xg, const double Ksac) const {
double a, b, c, tau, sq5, e1, e2, e3, e4, n1, n2, d1, d2, d3;

// Local copies of state and gating variables
const double V = X(0);
const double K_i = X(1);
const double Na_i = X(2);
const double Ca_i = X(3);
const double Ca_ss = X(4);
const double Ca_sr = X(5);
const double R_bar = X(6);
const double xr1 = Xg(0);
const double xr2 = Xg(1);
const double xs = Xg(2);
const double m = Xg(3);
const double h = Xg(4);
const double j = Xg(5);
const double d = Xg(6);
const double f = Xg(7);
const double f2 = Xg(8);
const double fcass = Xg(9);
const double s = Xg(10);
const double r = Xg(11);

const double RT = Rc * Tc / Fc;
const double E_K = RT * log(K_o / K_i);
const double E_Na = RT * log(Na_o / Na_i);
const double E_Ca = 0.5 * RT * log(Ca_o / Ca_i);
const double E_Ks = RT * log((K_o + p_KNa * Na_o) / (K_i + p_KNa * Na_i));

const double E_K_Ki = -RT / K_i;
const double E_Na_Nai = -RT / Na_i;
const double E_Ca_Cai = -RT / Ca_i / 2.0;
const double E_Ks_Ki = -RT / (K_i + p_KNa * Na_i);
const double E_Ks_Nai = p_KNa * E_Ks_Ki;

// I_Na: Fast sodium current
const double I_Na = G_Na * pow(m, 3.0) * h * j * (V - E_Na);
const double I_Na_V = G_Na * pow(m, 3.0) * h * j;
const double I_Na_Nai = I_Na_V * (-E_Na_Nai);

// I_to: transient outward current
const double I_to = G_to * r * s * (V - E_K);
const double I_to_V = G_to * r * s;
const double I_to_Ki = I_to_V * (-E_K_Ki);

// I_K1: inward rectifier outward current
e1 = exp(0.060 * (V - E_K - 200.0));
e2 = exp(2.E-40 * (V - E_K + 100.0));
e3 = exp(0.10 * (V - E_K - 10.0));
e4 = exp(-0.50 * (V - E_K));
a = 0.10 / (1.0 + e1);
b = (3.0 * e2 + e3) / (1.0 + e4);
tau = a / (a + b);
sq5 = sqrt(K_o / 5.40);
n1 = -6.E-30 * e1 / pow(1.0 + e1, 2.0);
n2 = (6.E-40 * e2 + 0.10 * e3 + 0.50 * b * e4) / (1.0 + e4);
n1 = (a + b) * n1 - a * (n1 + n2);
d1 = pow(a + b, 2.0);
const double I_K1 = G_K1 * sq5 * tau * (V - E_K);
const double I_K1_V = G_K1 * sq5 * (tau + (V - E_K) * n1 / d1);
const double I_K1_Ki = I_K1_V * (-E_K_Ki);

// I_Kr: rapid delayed rectifier current
const double I_Kr = G_Kr * sq5 * xr1 * xr2 * (V - E_K);
const double I_Kr_V = G_Kr * sq5 * xr1 * xr2;
const double I_Kr_Ki = I_Kr_V * (-E_K_Ki);

// I_Ks: slow delayed rectifier current
const double I_Ks = G_Ks * pow(xs, 2.0) * (V - E_Ks);
const double I_Ks_V = G_Ks * pow(xs, 2.0);
const double I_Ks_Ki = I_Ks_V * (-E_Ks_Ki);
const double I_Ks_Nai = I_Ks_V * (-E_Ks_Nai);

// I_CaL: L-type Ca current
a = 2.0 * (V - 15.0) / RT;
b = (0.250 * Ca_ss * exp(a) - Ca_o) / (exp(a) - 1.0);
c = G_CaL * d * f * f2 * fcass * (2.0 * a * Fc);
n1 = (exp(a) / RT) / (exp(a) - 1.0);
const double I_CaL = c * b;
const double I_CaL_V =
I_CaL / (V - 15.0) + n1 * (c * 0.50 * Ca_ss - 2.0 * I_CaL);
const double I_CaL_Cass = c * 0.250 * n1 * RT;

// I_NaCa: Na-Ca exchanger current
e1 = exp(gamma * V / RT);
e2 = exp((gamma - 1.0) * V / RT);
n1 = e1 * pow(Na_i, 3.0) * Ca_o - e2 * pow(Na_o, 3.0) * Ca_i * alpha;
d1 = pow(K_mNai, 3.0) + pow(Na_o, 3.0);
d2 = K_mCa + Ca_o;
d3 = 1.0 + K_sat * e2;
c = 1.0 / (d1 * d2 * d3);
const double I_NaCa = K_NaCa * n1 * c;

n1 = K_NaCa * c *
(e1 * pow(Na_i, 3.0) * Ca_o * (gamma / RT) -
e2 * pow(Na_o, 3.0) * Ca_i * alpha * ((gamma - 1.0) / RT));
n2 = I_NaCa * K_sat * ((gamma - 1.0) / RT) * e2 / d3;
const double I_NaCa_V = n1 - n2;
const double I_NaCa_Nai = K_NaCa * e1 * (3.0 * pow(Na_i, 2.0)) * Ca_o * c;
const double I_NaCa_Cai = -K_NaCa * e2 * pow(Na_o, 3.0) * alpha * c;

// I_NaK: Na-K pump current
e1 = exp(-0.10 * V / RT);
e2 = exp(-V / RT);
n1 = p_NaK * K_o * Na_i;
d1 = K_o + K_mK;
d2 = Na_i + K_mNa;
d3 = 1.0 + 0.12450 * e1 + 0.03530 * e2;
const double I_NaK = n1 / (d1 * d2 * d3);
n1 = (0.012450 * e1 + 0.03530 * e2) / RT;
const double I_NaK_V = I_NaK * n1 / d3;
const double I_NaK_Nai = I_NaK * K_mNa / (Na_i * d2);

// I_pCa: plateau Ca current
d1 = (K_pCa + Ca_i);
const double I_pCa = G_pCa * Ca_i / d1;
const double I_pCa_Cai = G_pCa * K_pCa / (d1 * d1);

// I_pK: plateau K current
e1 = exp((25.0 - V) / 5.980);
const double I_pK = G_pK * (V - E_K) / (1.0 + e1);
const double I_pK_V = (G_pK + I_pK * e1 / 5.980) / (1.0 + e1);
const double I_pK_Ki = G_pK * (-E_K_Ki) / (1.0 + e1);

// I_bCa: background Ca current
const double I_bCa = G_bCa * (V - E_Ca);
const double I_bCa_V = G_bCa;
const double I_bCa_Cai = G_bCa * (-E_Ca_Cai);

// I_bNa: background Na current
const double I_bNa = G_bNa * (V - E_Na);
const double I_bNa_V = G_bNa;
const double I_bNa_Nai = G_bNa * (-E_Na_Nai);

// I_leak: Sacroplasmic Reticulum Ca leak current
const double I_leak = V_leak * (Ca_sr - Ca_i);
const double I_leak_Cai = -V_leak;
const double I_leak_Casr = V_leak;

// I_up: Sacroplasmic Reticulum Ca pump current
d1 = 1.0 + pow(K_up / Ca_i, 2.0);
const double I_up = Vmax_up / d1;
const double I_up_Cai = (I_up / d1) * (2.0 * pow(K_up, 2.0) / pow(Ca_i, 3.0));

// I_rel: Ca induced Ca current (CICR)
n1 = max_sr - min_sr;
d1 = 1.0 + pow(EC / Ca_sr, 2.0);
const double k_casr = max_sr - (n1 / d1);
const double k1 = k1p / k_casr;
n2 = Ca_ss * 2.0;
d2 = k3 + k1 * n2;
const double O = k1 * R_bar * n2 / d2;
const double I_rel = V_rel * O * (Ca_sr - Ca_ss);

const double k_casr_sr =
(n1 / pow(d1, 2.0)) * (2.0 * pow(EC, 2.0) / pow(Ca_sr, 3.0));
// k_casr_sr = (n1 / (d1**2.0)) * (2.0*EC**2.0 / Ca_sr**3.0);
const double k1_casr = -k1p * k_casr_sr / pow(k_casr, 2.0);
const double O_Cass = 2.0 * k3 * O / (Ca_ss * d2);
const double O_Casr = k1_casr * n2 * (R_bar - O) / d2;
const double O_Rbar = k1 * n2 / d2;

const double I_rel_Cass = V_rel * (O_Cass * (Ca_sr - Ca_ss) - O);
const double I_rel_Casr = V_rel * (O_Casr * (Ca_sr - Ca_ss) + O);
const double I_rel_Rbar = V_rel * O_Rbar * (Ca_sr - Ca_ss);

// I_xfer: diffusive Ca current between Ca subspae and cytoplasm
const double I_xfer = V_xfer * (Ca_ss - Ca_i);
const double I_xfer_Cai = -V_xfer;
const double I_xfer_Cass = V_xfer;

// Compute Jacobian matrix
Array<double> Jac(X.size(), X.size());

c = Cm / (V_c * Fc);

// V
Jac(0, 0) = -(I_Na_V + I_to_V + I_K1_V + I_Kr_V + I_Ks_V + I_CaL_V +
I_NaCa_V + I_NaK_V + I_pK_V + I_bCa_V + I_bNa_V + Ksac);
Jac(0, 1) = -(I_to_Ki + I_K1_Ki + I_Kr_Ki + I_Ks_Ki + I_pK_Ki);
Jac(0, 2) = -(I_Na_Nai + I_Ks_Nai + I_NaCa_Nai + I_NaK_Nai + I_bNa_Nai);
Jac(0, 3) = -(I_NaCa_Cai + I_pCa_Cai + I_bCa_Cai);
Jac(0, 4) = -I_CaL_Cass;

// K_i
Jac(1, 0) = -c * (I_K1_V + I_to_V + I_Kr_V + I_Ks_V + I_pK_V - 2.0 * I_NaK_V);
Jac(1, 1) = -c * (I_K1_Ki + I_to_Ki + I_Kr_Ki + I_Ks_Ki + I_pK_Ki);
Jac(1, 2) = -c * (I_Ks_Nai - 2.0 * I_NaK_Nai);

// Na_i
Jac(2, 0) = -c * (I_Na_V + I_bNa_V + 3.0 * (I_NaK_V + I_NaCa_V));
Jac(2, 2) = -c * (I_Na_Nai + I_bNa_Nai + 3.0 * (I_NaK_Nai + I_NaCa_Nai));
Jac(2, 3) = -c * (3.0 * I_NaCa_Cai);

// Ca_i
n1 = (I_leak - I_up) * V_sr / V_c + I_xfer -
0.50 * c * (I_bCa + I_pCa - 2.0 * I_NaCa);
n2 = (I_leak_Cai - I_up_Cai) * V_sr / V_c + I_xfer_Cai -
0.50 * c * (I_bCa_Cai + I_pCa_Cai - 2.0 * I_NaCa_Cai);
d1 = 1.0 + K_bufc * Buf_c / pow(Ca_i + K_bufc, 2.0);
d2 = 2.0 * K_bufc * Buf_c / pow(Ca_i + K_bufc, 3.0);
Jac(3, 0) = -c * (I_bCa_V - 2.0 * I_NaCa_V) / 2.0 / d1;
Jac(3, 2) = c * I_NaCa_Nai / d1;
Jac(3, 3) = (n2 + n1 * d2 / d1) / d1;
Jac(3, 4) = I_xfer_Cass / d1;
Jac(3, 5) = (I_leak_Casr * V_sr / V_c) / d1;

// Ca_ss
a = Cm / (2.0 * Fc * V_ss);
b = V_sr / V_ss;
c = V_c / V_ss;
n1 = -a * I_CaL + b * I_rel - c * I_xfer;
n2 = -a * I_CaL_Cass + b * I_rel_Cass - c * I_xfer_Cass;
d1 = 1.0 + K_bufss * Buf_ss / pow(Ca_ss + K_bufss, 2.0);
d2 = 2.0 * K_bufss * Buf_ss / pow(Ca_ss + K_bufss, 3.0);
Jac(4, 0) = -a * I_CaL_V / d1;
Jac(4, 3) = -c * I_xfer_Cai / d1;
Jac(4, 4) = (n2 + n1 * d2 / d1) / d1;
Jac(4, 5) = b * I_rel_Casr / d1;
Jac(4, 6) = b * I_rel_Rbar / d1;

// Ca_sr
n1 = I_up - I_leak - I_rel;
n2 = -(I_leak_Casr + I_rel_Casr);
d1 = 1.0 + K_bufsr * Buf_sr / pow(Ca_sr + K_bufsr, 2.0);
d2 = 2.0 * K_bufsr * Buf_sr / pow(Ca_sr + K_bufsr, 3.0);
Jac(5, 3) = (I_up_Cai - I_leak_Cai) / d1;
Jac(5, 4) = -I_rel_Cass / d1;
Jac(5, 5) = (n2 + n1 * d2 / d1) / d1;
Jac(5, 6) = -I_rel_Rbar / d1;

// Rbar: ryanodine receptor
const double k2 = k2p * k_casr;
Jac(6, 4) = -k2 * R_bar;
Jac(6, 5) = -(k2p * k_casr_sr) * Ca_ss * R_bar;
Jac(6, 6) = -(k2 * Ca_ss + k4);

return Jac;
}

REGISTER_IONIC_MODEL("TTP", TTP);
5 changes: 0 additions & 5 deletions Code/Source/solver/ionic_ttp.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -326,11 +326,6 @@ class TTP : public IonicModel {
const Vector<double> &X, const Vector<double> &Xg,
const double I_stim,
const double I_sac) const override;

/// Model jacobian.
virtual Array<double> getj(const unsigned int zone_id,
const Vector<double> &X, const Vector<double> &Xg,
const double Ksac) const override;
};

#endif
5 changes: 0 additions & 5 deletions Code/Source/solver/read_files.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1119,11 +1119,6 @@ void read_cep_domain(Simulation* simulation, EquationParameters* eq_params, Doma
}
lDmn.cep.odes.tIntType = time_integration_type;

if ((lDmn.cep.odes.tIntType == TimeIntegrationType::CN2) &&
(lDmn.cep.cepType == ElectrophysiologyModelType::TTP)) {
throw std::runtime_error("[read_cep_domain] Implicit time integration for tenTusscher-Panfilov model can give unexpected results. Use FE or RK4 instead");
}

if (lDmn.cep.odes.tIntType == TimeIntegrationType::CN2) {
lDmn.cep.odes.maxItr = 5;
lDmn.cep.odes.absTol = 1e-8;
Expand Down
Loading