diff --git a/Code/Source/solver/ionic_ttp.cpp b/Code/Source/solver/ionic_ttp.cpp index efc6d6bea..c7cf58152 100644 --- a/Code/Source/solver/ionic_ttp.cpp +++ b/Code/Source/solver/ionic_ttp.cpp @@ -456,246 +456,5 @@ Vector TTP::getf(const unsigned int zone_id, const Vector &X, return dX; } -Array TTP::getj(const unsigned int zone_id, const Vector &X, - const Vector &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 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); \ No newline at end of file diff --git a/Code/Source/solver/ionic_ttp.h b/Code/Source/solver/ionic_ttp.h index afa5ad7b1..f23531d4c 100644 --- a/Code/Source/solver/ionic_ttp.h +++ b/Code/Source/solver/ionic_ttp.h @@ -326,11 +326,6 @@ class TTP : public IonicModel { const Vector &X, const Vector &Xg, const double I_stim, const double I_sac) const override; - - /// Model jacobian. - virtual Array getj(const unsigned int zone_id, - const Vector &X, const Vector &Xg, - const double Ksac) const override; }; #endif diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index fdfbe36ec..6c6723bce 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -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;