diff --git a/Code/Source/solver/CepModTtp.cpp b/Code/Source/solver/CepModTtp.cpp index c4a703660..be376ddaf 100644 --- a/Code/Source/solver/CepModTtp.cpp +++ b/Code/Source/solver/CepModTtp.cpp @@ -3,9 +3,35 @@ #include "CepModTtp.h" +#include "CmMod.h" +#include "Parameters.h" #include "mat_fun.h" #include +// Parameters discussed in https://doi.org/10.1152/ajpheart.00109.2006 +// Default parameters are for epicardium state (source: https://models.cellml.org/e/80d) +const TenTusscherPanfilovState TenTusscherPanfilovState::default_state = { + .V = -85.23, + .K_i = 136.89, + .Na_i = 8.6040, + .Ca_i = 1.26E-4, + .Ca_ss = 3.6E-4, + .Ca_sr = 3.64, + .R_bar = 0.9073, + .x_r1 = 6.21E-3, + .x_r2 = 0.4712, + .x_s = 9.5E-3, + .m = 1.72E-3, + .h = 0.7444, + .j = 0.7045, + .d = 3.373E-5, + .f = 0.7888, + .f2 = 0.9755, + .fcass = 0.9953, + .s = 0.999998, + .r = 2.42E-8 +}; + CepModTtp::CepModTtp() { } @@ -14,6 +40,7 @@ CepModTtp::~CepModTtp() { } + /// @brief Compute macroscopic fiber strain based on sacromere force-length relationship and calcium concentration void CepModTtp::actv_strn(const double c_Ca, const double I4f, const double dt, double& gf) { @@ -43,10 +70,8 @@ void CepModTtp::actv_strs(const double c_Ca, const double dt, double& Tact, doub /// @brief Compute currents and time derivatives of state variables /// -/// Note that is 'i' the myocardium zone id: 1, 2 or 3. -/// /// Reproduces Fortran 'GETF()'. -void CepModTtp::getf(const int i, const int nX, const int nG, const Vector& X, const Vector& Xg, Vector& dX, +void CepModTtp::getf(const int nX, const int nG, const Vector& X, const Vector& Xg, Vector& dX, const double I_stim, const double K_sac, Vector& RPAR) { // Local copies of state variables @@ -86,8 +111,7 @@ void CepModTtp::getf(const int i, const int nX, const int nG, const Vector& X, const Vector& Xg, +void CepModTtp::getj(const int nX, const int nG, const Vector& X, const Vector& Xg, Array& JAC, const double Ksac) { @@ -254,8 +278,8 @@ void CepModTtp::getj(const int i, const int nX, const int nG, const Vector& X, Vector& Xg ) +void CepModTtp::set_initial_conditions(const TTPInitialConditionsParameters& params) { - switch (imyo) { - - // epi - case 1: - // Initialize state variables - X(0) = -85.23; // V (units: mV) - X(1) = 136.89; // K_i (units: mM) - X(2) = 8.6040; // Na_i (units: mM) - X(3) = 1.26E-4; // Ca_i (units: mM) - X(4) = 3.6E-4; // Ca_ss (units: mM) - X(5) = 3.64; // Ca_sr (units: mM) - X(6) = 0.9073; // R' (dimensionless) - - // Initialize gating variables - Xg(0) = 6.21E-3; // x_r1 (dimensionless) - Xg(1) = 0.4712; // x_r2 (dimensionless) - Xg(2) = 9.5E-3; // x_s (dimensionless) - Xg(3) = 1.72E-3; // m (dimensionless) - Xg(4) = 0.7444; // h (dimensionless) - Xg(5) = 0.7045; // j (dimensionless) - Xg(6) = 3.373E-5; // d (dimensionless) - Xg(7) = 0.7888; // f (dimensionless) - Xg(8) = 0.9755; // f_2 (dimensionless) - Xg(9) = 0.9953; // f_cass (dimensionless) - Xg(10) = 0.999998; // s (dimensionless) - Xg(11) = 2.42E-8; // r (dimensionless) - break; - - // endo - case 2: - // Initialize state variables - X(0) = -86.709; // V (units: mV) - X(1) = 138.4; // K_i (units: mM) - X(2) = 10.355; // Na_i (units: mM) - X(3) = 1.3E-4; // Ca_i (units: mM) - X(4) = 3.6E-4; // Ca_ss (units: mM) - X(5) = 3.715; // Ca_sr (units: mM) - X(6) = 0.9068; // R' (dimensionless) - - // Initialize gating variables - Xg(0) = 4.48E-3; // x_r1 (dimensionless) - Xg(1) = 0.476; // x_r2 (dimensionless) - Xg(2) = 8.7E-3; // x_s (dimensionless) - Xg(3) = 1.55E-3; // m (dimensionless) - Xg(4) = 0.7573; // h (dimensionless) - Xg(5) = 0.7225; // j (dimensionless) - Xg(6) = 3.164E-5; // d (dimensionless) - Xg(7) = 0.8009; // f (dimensionless) - Xg(8) = 0.9778; // f_2 (dimensionless) - Xg(9) = 0.9953; // f_cass (dimensionless) - Xg(10) = 0.3212; // s (dimensionless) - Xg(11) = 2.235E-8; // r (dimensionless) - break; - - // mid-myo - case 3: - // Initialize state variables - X(0) = -85.423; // V (units: mV) - X(1) = 138.52; // K_i (units: mM) - X(2) = 10.132; // Na_i (units: mM) - X(3) = 1.53E-4; // Ca_i (units: mM) - X(4) = 4.2E-4; // Ca_ss (units: mM) - X(5) = 4.272; // Ca_sr (units: mM) - X(6) = 0.8978; // R' (dimensionless) - - // Initialize gating variables - Xg(0) = 1.65E-2; // x_r1 (dimensionless) - Xg(1) = 0.4730; // x_r2 (dimensionless) - Xg(2) = 1.74E-2; // x_s (dimensionless) - Xg(3) = 1.65E-3; // m (dimensionless) - Xg(4) = 0.7490; // h (dimensionless) - Xg(5) = 0.6788; // j (dimensionless) - Xg(6) = 3.288E-5; // d (dimensionless) - Xg(7) = 0.7026; // f (dimensionless) - Xg(8) = 0.9526; // f_2 (dimensionless) - Xg(9) = 0.9942; // f_cass (dimensionless) - Xg(10) = 0.999998; // s (dimensionless) - Xg(11) = 2.347E-8; // r (dimensionless) - break; - } - + auto& s = params.initial_states; + initial_state.V = s.V.value(); + initial_state.K_i = s.K_i.value(); + initial_state.Na_i = s.Na_i.value(); + initial_state.Ca_i = s.Ca_i.value(); + initial_state.Ca_ss = s.Ca_ss.value(); + initial_state.Ca_sr = s.Ca_sr.value(); + initial_state.R_bar = s.R_bar.value(); + + auto& g = params.gating_variables; + initial_state.x_r1 = g.x_r1_rectifier.value(); + initial_state.x_r2 = g.x_r2_rectifier.value(); + initial_state.x_s = g.x_s_rectifier.value(); + initial_state.m = g.m_fast_Na.value(); + initial_state.h = g.h_fast_Na.value(); + initial_state.j = g.j_fast_Na.value(); + initial_state.d = g.d_slow_in.value(); + initial_state.f = g.f_slow_in.value(); + initial_state.f2 = g.f2_slow_in.value(); + initial_state.fcass = g.fcass_slow_in.value(); + initial_state.s = g.s_out.value(); + initial_state.r = g.r_out.value(); + + user_initial_state = true; } -void CepModTtp::init(const int imyo, const int nX, const int nG, Vector& X, Vector& Xg, - Vector& X0, Vector& Xg0) +void CepModTtp::distribute_conductance(const CmMod& cm_mod, const cmType& cm) { - init(imyo, nX, nG, X, Xg); + cm.bcast(cm_mod, &G_Na); + cm.bcast(cm_mod, &G_CaL); + cm.bcast(cm_mod, &G_Kr); + cm.bcast(cm_mod, &G_Ks); + cm.bcast(cm_mod, &G_to); +} - if (X0.size() != 0) { - X = X0; +void CepModTtp::distribute_initial_state(const CmMod& cm_mod, const cmType& cm) +{ + cm.bcast(cm_mod, &user_initial_state); + + if (user_initial_state) { + cm.bcast(cm_mod, &initial_state.V); + cm.bcast(cm_mod, &initial_state.K_i); + cm.bcast(cm_mod, &initial_state.Na_i); + cm.bcast(cm_mod, &initial_state.Ca_i); + cm.bcast(cm_mod, &initial_state.Ca_ss); + cm.bcast(cm_mod, &initial_state.Ca_sr); + cm.bcast(cm_mod, &initial_state.R_bar); + cm.bcast(cm_mod, &initial_state.x_r1); + cm.bcast(cm_mod, &initial_state.x_r2); + cm.bcast(cm_mod, &initial_state.x_s); + cm.bcast(cm_mod, &initial_state.m); + cm.bcast(cm_mod, &initial_state.h); + cm.bcast(cm_mod, &initial_state.j); + cm.bcast(cm_mod, &initial_state.d); + cm.bcast(cm_mod, &initial_state.f); + cm.bcast(cm_mod, &initial_state.f2); + cm.bcast(cm_mod, &initial_state.fcass); + cm.bcast(cm_mod, &initial_state.s); + cm.bcast(cm_mod, &initial_state.r); } +} - if (Xg0.size() != 0) { - Xg = Xg0; +void CepModTtp::init(const int nX, const int nG, Vector& X, Vector& Xg) +{ + if (!user_initial_state) { + initial_state = TenTusscherPanfilovState::default_state; } + copy_state_to_vectors(X, Xg); } + /// @brief Time integration performed using Crank-Nicholson method void CepModTtp::integ_cn2(const int imyo, const int nX, const int nG, Vector& Xn, Vector& Xg, const double Ts, const double dt, const double Istim, const double Ksac, Vector& IPAR, Vector& RPAR) @@ -553,7 +551,7 @@ void CepModTtp::integ_cn2(const int imyo, const int nX, const int nG, Vector fn(nX); - getf(imyo, nX, nG, Xn, Xg, fn, Istim, Ksac, RPAR); + getf(nX, nG, Xn, Xg, fn, Istim, Ksac, RPAR); int k = 0; auto Xk = Xn; @@ -566,7 +564,7 @@ void CepModTtp::integ_cn2(const int imyo, const int nX, const int nG, Vector fk(nX); - getf(imyo, nX, nG, Xn, Xg, fn, Istim, Ksac, RPAR); + getf(nX, nG, Xn, Xg, fn, Istim, Ksac, RPAR); auto rK = Xk - Xn - 0.5*dt*(fk + fn); @@ -589,7 +587,7 @@ void CepModTtp::integ_cn2(const int imyo, const int nX, const int nG, Vector JAC(nX,nX); - getj(imyo, nX, nG, Xk, Xg, JAC, Ksac); + getj(nX, nG, Xk, Xg, JAC, Ksac); JAC = Im - 0.5*dt*JAC; JAC = mat_fun::mat_inv(JAC, nX); @@ -600,7 +598,7 @@ void CepModTtp::integ_cn2(const int imyo, const int nX, const int nG, Vector f(nX); // Get time derivatives (RHS) - getf(imyo, nX, nG, X, Xg, f, Istim, Ksac, RPAR); + getf(nX, nG, X, Xg, f, Istim, Ksac, RPAR); //CALL TTP_GETF(imyo, nX, nG, X, Xg, f, Istim, Ksac, RPAR) // Update gating variables @@ -633,7 +631,7 @@ void CepModTtp::integ_rk(const int imyo, const int nX, const int nG, Vector& X, Vector& Xg) +void CepModTtp::update_g(const int imyo, const double dt, const int n, const int nG, const Vector& X, Vector& Xg) { V = X(0); Ca_ss = X(4); @@ -774,10 +772,10 @@ void CepModTtp::update_g(const int i, const double dt, const int n, const int nG Xg(9) = fcassi - (fcassi - fcass)*exp(-dt/tau); // s: inactivation gate for I_to - if (i == 1 || i == 3) { + if (imyo == 1 || imyo == 3) { si = 1.0/(1.0 + exp((20.0+V)/5.0)); tau = 85.0*exp(-pow(V+45.0,2.0) / 320.0) + 5.0/(1.0+exp((V-20.0)/5.0)) + 3.0; - } else if (i == 2) { + } else if (imyo == 2) { si = 1.0/(1.0 + exp((28.0+V)/5.0)); tau = 1000.0*exp(-pow(V+67.0,2.0) /1000.0) + 8.0; //tau = 1000.0*exp(-((V+67.0)**2.0) /1000.0) + 8.0; @@ -791,5 +789,26 @@ void CepModTtp::update_g(const int i, const double dt, const int n, const int nG Xg(11) = ri - (ri - r)*exp(-dt/tau); } - - +void CepModTtp::copy_state_to_vectors(Vector& X, Vector& Xg) const +{ + X(0) = initial_state.V; + X(1) = initial_state.K_i; + X(2) = initial_state.Na_i; + X(3) = initial_state.Ca_i; + X(4) = initial_state.Ca_ss; + X(5) = initial_state.Ca_sr; + X(6) = initial_state.R_bar; + + Xg(0) = initial_state.x_r1; + Xg(1) = initial_state.x_r2; + Xg(2) = initial_state.x_s; + Xg(3) = initial_state.m; + Xg(4) = initial_state.h; + Xg(5) = initial_state.j; + Xg(6) = initial_state.d; + Xg(7) = initial_state.f; + Xg(8) = initial_state.f2; + Xg(9) = initial_state.fcass; + Xg(10) = initial_state.s; + Xg(11) = initial_state.r; +} diff --git a/Code/Source/solver/CepModTtp.h b/Code/Source/solver/CepModTtp.h index 181a7b7f0..d8d5b9dc8 100644 --- a/Code/Source/solver/CepModTtp.h +++ b/Code/Source/solver/CepModTtp.h @@ -8,6 +8,10 @@ #include "Vector.h" #include +class CmMod; +class cmType; +class TTPInitialConditionsParameters; + #include #include @@ -17,10 +21,40 @@ T& make_ref(T&& x) { return x; } /// @brief This module defines data structures for ten Tusscher-Panfilov /// epicardial cellular activation model for cardiac electrophysiology /// -/// The classes defined here duplicate the data structures in the Fortran TPPMOD module defined -/// in CEPMOD_TTP.f and PARAMS_TPP.f files. +/// The classes defined here duplicate the data structures in the Fortran TPPMOD module defined +/// in CEPMOD_TTP.f and PARAMS_TPP.f files. + +/// @brief Stores all TTP state and gating variables in a single typed struct. +class TenTusscherPanfilovState { + public: + double V; + double K_i; + double Na_i; + double Ca_i; + double Ca_ss; + double Ca_sr; + double R_bar; + double x_r1; + double x_r2; + double x_s; + double m; + double h; + double j; + double d; + double f; + double f2; + double fcass; + double s; + double r; + + static const TenTusscherPanfilovState default_state; +}; + class CepModTtp { + private: + TenTusscherPanfilovState initial_state; + public: CepModTtp(); ~CepModTtp(); @@ -29,7 +63,7 @@ class CepModTtp //-------------------------------------------------------------------- // // Constants for TenTusscher-Panfilov Ventricular Myocyte Model. -// +// Default values are for the epicardium state (source: https://models.cellml.org/e/80d) //-------------------------------------------------------------------- // Default model parameters @@ -75,20 +109,14 @@ class CepModTtp /// Maximal I_K1 conductance [nS/pF] double G_K1 = 5.405; - /// Maximal epicardial I_to conductance [nS/pF] - Vector G_to = {0.294, 0.073, 0.294}; + /// Maximal I_to conductance [nS/pF] + double G_to = 0.294; /// Maximal I_Kr conductance [nS/pF] double G_Kr = 0.153; -// G_Kr for spiral wave breakup -// double G_Kr = 0.172; // units: nS/pF - - /// Maximal epicardial I_Ks conductance [nS/pF] - Vector G_Ks = {0.392, 0.392, 0.098}; - -// G_Ks for spiral wave breakup (epi) -// double G_Ks(3) = (/0.441, 0.392_RKIND, 0.098_RKIND/) + /// Maximal I_Ks conductance [nS/pF] + double G_Ks = 0.392; /// Relative I_Ks permeability to Na [-] double p_KNa = 3.E-2; @@ -126,15 +154,9 @@ class CepModTtp /// Maximal I_pK conductance [nS/pF] double G_pK = 1.46E-2; -// G_pK for spiral wave breakup -// double G_pK = 2.19E-3; // units: nS/pF - /// Maximal I_pCa conductance [pA/pF] double G_pCa = 0.1238; -// G_pCa for spiral wave breakup -// double G_pCa = 0.8666; // units: pA/pF - /// Half-saturation constant of I_pCa [mM] double K_pCa = 5.E-4; @@ -367,19 +389,21 @@ class CepModTtp double I_xfer_Cai, I_xfer_Cass; double k_casr_sr, k1_casr, O_Casr, O_Cass, O_Rbar; +// Flag for user defined initial conditions + bool user_initial_state = false; + void actv_strn(const double c_Ca, const double I4f, const double dt, double& gf); void actv_strs(const double c_Ca, const double dt, double& Tact, double& epsX); - void getf(const int i, const int nX, const int nG, const Vector& X, const Vector& Xg, + void getf(const int nX, const int nG, const Vector& X, const Vector& Xg, Vector& dX, const double I_stim, const double K_sac, Vector& RPAR); - void getj(const int i, const int nX, const int nG, const Vector& X, const Vector& Xg, + void getj(const int nX, const int nG, const Vector& X, const Vector& Xg, Array& JAC, const double Ksac); - void init(const int imyo, const int nX, const int nG, Vector& X, Vector& Xg); + void init(const int nX, const int nG, Vector& X, Vector& Xg); - void init(const int imyo, const int nX, const int nG, Vector& X, Vector& Xg, - Vector& X0, Vector& Xg0); + void set_initial_conditions(const TTPInitialConditionsParameters& params); void integ_cn2(const int imyo, const int nX, const int nG, Vector& X, Vector& Xg, const double Ts, const double dt, const double Istim, const double Ksac, @@ -391,9 +415,15 @@ class CepModTtp void integ_rk(const int imyo, const int nX, const int nG, Vector& X, Vector& Xg, const double Ts, const double dt, const double Istim, const double Ksac, Vector& RPAR); - void update_g(const int i, const double dt, const int n, const int nG, const Vector& X, + void update_g(const int imyo, const double dt, const int n, const int nG, const Vector& X, Vector& Xg); + void copy_state_to_vectors(Vector& X, Vector& Xg) const; + + void distribute_conductance(const CmMod& cm_mod, const cmType& cm); + + void distribute_initial_state(const CmMod& cm_mod, const cmType& cm); + }; #endif diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index d308e89b7..e5dc908c0 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -1666,9 +1666,12 @@ void DomainParameters::print_parameters() stimulus.print_parameters(); + ttp_initial_conditions.print_parameters(); + fluid_viscosity.print_parameters(); solid_viscosity.print_parameters(); + } //------------ @@ -1711,6 +1714,9 @@ void DomainParameters::set_values(tinyxml2::XMLElement* domain_elem, bool from_e } else if (name == StimulusParameters::xml_element_name_) { stimulus.set_values(item); + } else if (name == TTPInitialConditionsParameters::xml_element_name_) { + ttp_initial_conditions.set_values(item); + } else if (name == FluidViscosityParameters::xml_element_name_ || name == SolidViscosityParameters::xml_element_name_) { auto eq_type = consts::equation_name_to_type.at(equation.value()); if (eq_type == consts::EquationType::phys_fluid || eq_type == consts::EquationType::phys_CMM || eq_type == consts::EquationType::phys_stokes) { @@ -1758,6 +1764,187 @@ void DomainParameters::set_values(tinyxml2::XMLElement* domain_elem, bool from_e */ } +////////////////////////////////////////////////////////// +// TTPInitialConditionsParameters // +////////////////////////////////////////////////////////// + +const std::string TTPInitialConditionsParameters::xml_element_name_ = "TTP_initial_conditions"; + +TTPInitialConditionsParameters::TTPInitialConditionsParameters() +{ +} + +void TTPInitialConditionsParameters::print_parameters() +{ + if (value_set) { + std::cout << std::endl; + std::cout << "TTP Initial Conditions Parameters" << std::endl; + std::cout << "---------------------------------" << std::endl; + initial_states.print_parameters(); + gating_variables.print_parameters(); + } +} + +void TTPInitialConditionsParameters::set_values(tinyxml2::XMLElement* xml_elem) +{ + using namespace tinyxml2; + std::string error_msg = "Unknown " + xml_element_name_ + " XML element '"; + + auto item = xml_elem->FirstChildElement(); + + while (item != nullptr) { + auto name = std::string(item->Value()); + + if (name == TTPInitialStatesParameters::xml_element_name_) { + initial_states.set_values(item); + value_set = true; + + } else if (name == TTPGatingVariablesParameters::xml_element_name_) { + gating_variables.set_values(item); + value_set = true; + + } else { + throw std::runtime_error(error_msg + name + "'."); + } + + item = item->NextSiblingElement(); + } + + if (!initial_states.defined()) { + throw std::runtime_error(xml_element_name_ + " requires an '" + + TTPInitialStatesParameters::xml_element_name_ + "' XML section."); + } + + if (!gating_variables.defined()) { + throw std::runtime_error(xml_element_name_ + " requires a '" + + TTPGatingVariablesParameters::xml_element_name_ + "' XML section."); + } +} + +////////////////////////////////////////////////////////// +// TTPInitialStatesParameters // +////////////////////////////////////////////////////////// + +const std::string TTPInitialStatesParameters::xml_element_name_ = "Initial_states"; + +TTPInitialStatesParameters::TTPInitialStatesParameters() +{ + bool required = true; + + set_parameter("V", -85.23, required, V); + set_parameter("K_i", 136.89, required, K_i); + set_parameter("Na_i", 8.6040, required, Na_i); + set_parameter("Ca_i", 1.26E-4, required, Ca_i); + set_parameter("Ca_ss", 3.6E-4, required, Ca_ss); + set_parameter("Ca_sr", 3.64, required, Ca_sr); + set_parameter("R_bar", 0.9073, required, R_bar); +} + +void TTPInitialStatesParameters::print_parameters() +{ + if (value_set) { + std::cout << " Initial States:" << std::endl; + auto params_name_value = get_parameter_list(); + for (auto& [key, value] : params_name_value) { + std::cout << " " << key << ": " << value << std::endl; + } + } +} + +void TTPInitialStatesParameters::set_values(tinyxml2::XMLElement* xml_elem) +{ + using namespace tinyxml2; + std::string error_msg = "Unknown " + xml_element_name_ + " XML element '"; + + auto item = xml_elem->FirstChildElement(); + + while (item != nullptr) { + auto name = std::string(item->Value()); + + if (item->GetText() != nullptr) { + auto value = item->GetText(); + try { + set_parameter_value(name, value); + value_set = true; + } catch (const std::bad_function_call& exception) { + throw std::runtime_error(error_msg + name + "'."); + } + } else { + throw std::runtime_error(error_msg + name + "'."); + } + + item = item->NextSiblingElement(); + } + + check_required(); +} + +////////////////////////////////////////////////////////// +// TTPGatingVariablesParameters // +////////////////////////////////////////////////////////// + +const std::string TTPGatingVariablesParameters::xml_element_name_ = "Gating_variables"; + +TTPGatingVariablesParameters::TTPGatingVariablesParameters() +{ + bool required = true; + + set_parameter("x_r1_rectifier", 6.21E-3, required, x_r1_rectifier); + set_parameter("x_r2_rectifier", 0.4712, required, x_r2_rectifier); + set_parameter("x_s_rectifier", 9.5E-3, required, x_s_rectifier); + + set_parameter("m_fast_Na", 1.72E-3, required, m_fast_Na); + set_parameter("h_fast_Na", 0.7444, required, h_fast_Na); + set_parameter("j_fast_Na", 0.7045, required, j_fast_Na); + + set_parameter("d_slow_in", 3.373E-5, required, d_slow_in); + set_parameter("f_slow_in", 0.7888, required, f_slow_in); + set_parameter("f2_slow_in", 0.9755, required, f2_slow_in); + set_parameter("fcass_slow_in", 0.9953, required, fcass_slow_in); + + set_parameter("s_out", 0.999998, required, s_out); + set_parameter("r_out", 2.42E-8, required, r_out); +} + +void TTPGatingVariablesParameters::print_parameters() +{ + if (value_set) { + std::cout << " Gating Variables:" << std::endl; + auto params_name_value = get_parameter_list(); + for (auto& [key, value] : params_name_value) { + std::cout << " " << key << ": " << value << std::endl; + } + } +} + +void TTPGatingVariablesParameters::set_values(tinyxml2::XMLElement* xml_elem) +{ + using namespace tinyxml2; + std::string error_msg = "Unknown " + xml_element_name_ + " XML element '"; + + auto item = xml_elem->FirstChildElement(); + + while (item != nullptr) { + auto name = std::string(item->Value()); + + if (item->GetText() != nullptr) { + auto value = item->GetText(); + try { + set_parameter_value(name, value); + value_set = true; + } catch (const std::bad_function_call& exception) { + throw std::runtime_error(error_msg + name + "'."); + } + } else { + throw std::runtime_error(error_msg + name + "'."); + } + + item = item->NextSiblingElement(); + } + + check_required(); +} + ////////////////////////////////////////////////////////// // DirectionalDistributionParameters // ////////////////////////////////////////////////////////// @@ -3137,4 +3324,3 @@ void LinearSolverParameters::set_values(tinyxml2::XMLElement* xml_elem) } } - diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index d4db1580b..45ba95dcb 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1180,6 +1180,81 @@ class FiberReinforcementStressParameters : public ParameterLists bool value_set = false; }; +/// @brief Stores parameters for the 'Gating_variables' XML element +/// under TTP_initial_conditions. +class TTPGatingVariablesParameters : public ParameterLists +{ + public: + TTPGatingVariablesParameters(); + + static const std::string xml_element_name_; + + bool defined() const { return value_set; }; + void print_parameters(); + void set_values(tinyxml2::XMLElement* xml_elem); + + Parameter x_r1_rectifier; + Parameter x_r2_rectifier; + Parameter x_s_rectifier; + + Parameter m_fast_Na; + Parameter h_fast_Na; + Parameter j_fast_Na; + + Parameter d_slow_in; + Parameter f_slow_in; + Parameter f2_slow_in; + Parameter fcass_slow_in; + + Parameter s_out; + Parameter r_out; + + bool value_set = false; +}; + +/// @brief Stores parameters for the 'Initial_states' XML element +/// under TTP_initial_conditions. +class TTPInitialStatesParameters : public ParameterLists +{ + public: + TTPInitialStatesParameters(); + + static const std::string xml_element_name_; + + bool defined() const { return value_set; }; + void print_parameters(); + void set_values(tinyxml2::XMLElement* xml_elem); + + Parameter V; + Parameter K_i; + Parameter Na_i; + Parameter Ca_i; + Parameter Ca_ss; + Parameter Ca_sr; + Parameter R_bar; + + bool value_set = false; +}; + +/// @brief Stores parameters for the 'TTP_initial_conditions' XML element +/// under Domain. +class TTPInitialConditionsParameters : public ParameterLists +{ + public: + TTPInitialConditionsParameters(); + + static const std::string xml_element_name_; + + bool defined() const { return value_set; }; + void print_parameters(); + void set_values(tinyxml2::XMLElement* xml_elem); + + TTPInitialStatesParameters initial_states; + TTPGatingVariablesParameters gating_variables; + + bool value_set = false; +}; + /// @brief The DomainParameters class stores parameters for the XML /// 'Domain' element to specify properties for solving equations. /// @@ -1209,6 +1284,7 @@ class DomainParameters : public ParameterLists StimulusParameters stimulus; FluidViscosityParameters fluid_viscosity; SolidViscosityParameters solid_viscosity; + TTPInitialConditionsParameters ttp_initial_conditions; // Attributes. Parameter id; diff --git a/Code/Source/solver/cep_ion.cpp b/Code/Source/solver/cep_ion.cpp index 2cf78b836..4355a4960 100644 --- a/Code/Source/solver/cep_ion.cpp +++ b/Code/Source/solver/cep_ion.cpp @@ -42,7 +42,7 @@ void cep_init(Simulation* simulation) } if (com_mod.dmnId.size() != 0) { - Vector sA(tnNo); + Vector sA(tnNo); Array sF(nXion,tnNo); for (int a = 0; a < tnNo; a++) { @@ -131,7 +131,7 @@ void cep_init_l(cepModelType& cep, int nX, int nG, Vector& X, Vector*,double*> simple_ttp_params{ - {&domain_params->G_Na, &lDmn.cep.ttp.G_Na}, - {&domain_params->G_Kr, &lDmn.cep.ttp.G_Kr}, - {&domain_params->G_CaL, &lDmn.cep.ttp.G_CaL} + {&domain_params->G_Na, &lDmn.cep.ttp.G_Na}, + {&domain_params->G_Kr, &lDmn.cep.ttp.G_Kr}, + {&domain_params->G_CaL, &lDmn.cep.ttp.G_CaL}, + {&domain_params->G_Ks, &lDmn.cep.ttp.G_Ks}, + {&domain_params->G_to, &lDmn.cep.ttp.G_to}, }; for (auto& [param, value] : simple_ttp_params) { @@ -1013,12 +1014,8 @@ void read_cep_domain(Simulation* simulation, EquationParameters* eq_params, Doma } } - // Array params - 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(); + if (model_type == ElectrophysiologyModelType::TTP) { + lDmn.cep.ttp.set_initial_conditions(domain_params->ttp_initial_conditions); } // Set stimulus parameters. diff --git a/tests/cases/cep/cable_TTP_1d/solver.xml b/tests/cases/cep/cable_TTP_1d/solver.xml index 9441f5035..bacef4b3b 100644 --- a/tests/cases/cep/cable_TTP_1d/solver.xml +++ b/tests/cases/cep/cable_TTP_1d/solver.xml @@ -36,12 +36,14 @@ TTP + ../ttp_parameters/ttp_epicardium_parameters.xml 0.15432 RK TTP + ../ttp_parameters/ttp_epicardium_parameters.xml 0.15432 RK diff --git a/tests/cases/cep/cylinder_purkinje_1d3d/solver.xml b/tests/cases/cep/cylinder_purkinje_1d3d/solver.xml index 35844bcb0..cef558fd5 100644 --- a/tests/cases/cep/cylinder_purkinje_1d3d/solver.xml +++ b/tests/cases/cep/cylinder_purkinje_1d3d/solver.xml @@ -61,7 +61,7 @@ TTP - myocardium + ../ttp_parameters/ttp_midmyocardium_parameters.xml 0.05 0.05 RK @@ -69,14 +69,14 @@ TTP - pfib + ../ttp_parameters/ttp_endocardium_parameters.xml 1.1 RK TTP - pfib + ../ttp_parameters/ttp_endocardium_parameters.xml 1.1 RK diff --git a/tests/cases/cep/niederer_benchmark_ECGs_quadrature/solver_GMRES_FE_epicardium_TTP.xml b/tests/cases/cep/niederer_benchmark_ECGs_quadrature/solver_GMRES_FE_epicardium_TTP.xml index 6afd15d44..09b75e5aa 100644 --- a/tests/cases/cep/niederer_benchmark_ECGs_quadrature/solver_GMRES_FE_epicardium_TTP.xml +++ b/tests/cases/cep/niederer_benchmark_ECGs_quadrature/solver_GMRES_FE_epicardium_TTP.xml @@ -76,13 +76,11 @@ 14.838 3.98E-5 0.153 - 0.392 - 0.294 0.012571 0.082715 - epicardium + ../ttp_parameters/ttp_epicardium_parameters.xml FE @@ -93,13 +91,11 @@ 14.838 3.98E-5 0.153 - 0.392 - 0.294 0.012571 0.082715 - epicardium + ../ttp_parameters/ttp_epicardium_parameters.xml -35.714 diff --git a/tests/cases/cep/niederer_benchmark_ECGs_quadrature/solver_GMRES_FE_pfib_AP.xml b/tests/cases/cep/niederer_benchmark_ECGs_quadrature/solver_GMRES_FE_pfib_AP.xml index acd86172b..f782bfcc4 100644 --- a/tests/cases/cep/niederer_benchmark_ECGs_quadrature/solver_GMRES_FE_pfib_AP.xml +++ b/tests/cases/cep/niederer_benchmark_ECGs_quadrature/solver_GMRES_FE_pfib_AP.xml @@ -76,8 +76,6 @@ 0.012571 0.082715 - pfib - FE @@ -87,8 +85,6 @@ 0.012571 0.082715 - pfib - -35.714 0.0 diff --git a/tests/cases/cep/purkinje/solver.xml b/tests/cases/cep/purkinje/solver.xml index 395aee00b..0f26e1e88 100644 --- a/tests/cases/cep/purkinje/solver.xml +++ b/tests/cases/cep/purkinje/solver.xml @@ -44,14 +44,14 @@ TTP - pfib + ../ttp_parameters/ttp_endocardium_parameters.xml 1.1 RK TTP - pfib + ../ttp_parameters/ttp_endocardium_parameters.xml 1.1 RK diff --git a/tests/cases/cep/slab_domains/solver_dat.xml b/tests/cases/cep/slab_domains/solver_dat.xml index fd626d280..8f138d899 100644 --- a/tests/cases/cep/slab_domains/solver_dat.xml +++ b/tests/cases/cep/slab_domains/solver_dat.xml @@ -50,7 +50,7 @@ 29.68 7.96e-5 0.31 - epi + ../ttp_parameters/ttp_epicardium_parameters.xml 0.033192 RK @@ -60,7 +60,7 @@ 7.42 1.99e-5 0.08 - epi + ../ttp_parameters/ttp_epicardium_parameters.xml 0.033192 RK @@ -72,8 +72,8 @@ 7.96e-5 0.31 - epi - + ../ttp_parameters/ttp_epicardium_parameters.xml + 0.033192 diff --git a/tests/cases/cep/slab_domains/solver_vtu.xml b/tests/cases/cep/slab_domains/solver_vtu.xml index e82c3e8b2..2f9c3adf0 100644 --- a/tests/cases/cep/slab_domains/solver_vtu.xml +++ b/tests/cases/cep/slab_domains/solver_vtu.xml @@ -50,7 +50,7 @@ 29.68 7.96e-5 0.31 - epi + ../ttp_parameters/ttp_epicardium_parameters.xml 0.033192 RK @@ -60,7 +60,7 @@ 7.42 1.99e-5 0.08 - epi + ../ttp_parameters/ttp_epicardium_parameters.xml 0.033192 RK @@ -72,7 +72,7 @@ 7.96e-5 0.31 - epi + ../ttp_parameters/ttp_epicardium_parameters.xml 0.033192 diff --git a/tests/cases/cep/spiral_TTP_2d/README.md b/tests/cases/cep/spiral_TTP_2d/README.md new file mode 100644 index 000000000..08231c6c4 --- /dev/null +++ b/tests/cases/cep/spiral_TTP_2d/README.md @@ -0,0 +1,9 @@ +# Spiral Wave for TTP Model + +Another approach to spiral wave initialization is to 'rig' the simulation domain such that there are 4 sub-domains which are all initialized to different sets of ionic and gating variable states. A visualization of this simulation case is shown below + +

+ +

+ +Each set of initial conditions parameters is defined in a seperate xml file that is included in the solver.xml via the `` parameter. diff --git a/tests/cases/cep/spiral_TTP_2d/SpiralWave.jpg b/tests/cases/cep/spiral_TTP_2d/SpiralWave.jpg new file mode 100644 index 000000000..a9e4aaf88 Binary files /dev/null and b/tests/cases/cep/spiral_TTP_2d/SpiralWave.jpg differ diff --git a/tests/cases/cep/spiral_TTP_2d/domain-1_EPmodel.xml b/tests/cases/cep/spiral_TTP_2d/domain-1_EPmodel.xml new file mode 100644 index 000000000..db64441ae --- /dev/null +++ b/tests/cases/cep/spiral_TTP_2d/domain-1_EPmodel.xml @@ -0,0 +1,41 @@ + + + + 14.838 + 3.98E-5 + 0.153 + 0.392 + 0.294 + + 0.1 + + epicardium + + + + 18.167 + 136.897 + 8.6105 + 1.2592E-4 + 3.6181E-4 + 3.6399 + 0.9078 + + + + 6.368E-3 + 0.36774 + 9.77E-3 + 0.9338 + 0.5276 + 0.6837 + 1.9520E-2 + 0.7896 + 0.9752 + 0.9954 + 0.9941 + 4.2757E-5 + + + + \ No newline at end of file diff --git a/tests/cases/cep/spiral_TTP_2d/domain-2_EPmodel.xml b/tests/cases/cep/spiral_TTP_2d/domain-2_EPmodel.xml new file mode 100644 index 000000000..9ca684de3 --- /dev/null +++ b/tests/cases/cep/spiral_TTP_2d/domain-2_EPmodel.xml @@ -0,0 +1,41 @@ + + + + 14.838 + 3.98E-5 + 0.153 + 0.392 + 0.294 + + 0.1 + + epicardium + + + + -85.423 + 138.52 + 10.132 + 1.53E-4 + 4.2E-4 + 4.272 + 0.8978 + + + + 1.65E-2 + 0.4730 + 1.74E-2 + 1.65E-3 + 0.7490 + 0.6788 + 3.288E-5 + 0.7026 + 0.9526 + 0.9942 + 0.999998 + 2.347E-8 + + + + \ No newline at end of file diff --git a/tests/cases/cep/spiral_TTP_2d/domain-3_EPmodel.xml b/tests/cases/cep/spiral_TTP_2d/domain-3_EPmodel.xml new file mode 100644 index 000000000..49b0a8626 --- /dev/null +++ b/tests/cases/cep/spiral_TTP_2d/domain-3_EPmodel.xml @@ -0,0 +1,41 @@ + + + + 14.838 + 3.98E-5 + 0.153 + 0.392 + 0.294 + + 0.1 + + epicardium + + + + -44.756 + 136.88 + 8.5907 + 3.0084E-4 + 2.5462E-3 + 3.5292 + 0.7235 + + + + 0.9739 + 0.1239 + 0.1766 + 0.6379 + 2.3557E-4 + 3.147E-5 + 9.7186E-3 + 0.1754 + 0.3475 + 0.9296 + 0.4095 + 7.857E-4 + + + + \ No newline at end of file diff --git a/tests/cases/cep/spiral_TTP_2d/domain-4_EPmodel.xml b/tests/cases/cep/spiral_TTP_2d/domain-4_EPmodel.xml new file mode 100644 index 000000000..04f1ebd4b --- /dev/null +++ b/tests/cases/cep/spiral_TTP_2d/domain-4_EPmodel.xml @@ -0,0 +1,41 @@ + + + + 14.838 + 3.98E-5 + 0.153 + 0.392 + 0.294 + + 0.1 + + epicardium + + + + -11.612 + 136.88 + 8.5916 + 3.486E-4 + 1.312E-2 + 3.4619 + 0.6875 + + + + 9.882E-1 + 3.707E-2 + 1.7887E-1 + 9.8692E-1 + 9.184E-8 + 6.291E-8 + 0.4147 + 0.15883 + 0.3341 + 0.9036 + 7.9355E-2 + 1.0429E-2 + + + + \ No newline at end of file diff --git a/tests/cases/cep/spiral_TTP_2d/mesh/domain_ids.dat b/tests/cases/cep/spiral_TTP_2d/mesh/domain_ids.dat new file mode 100644 index 000000000..bf695480e --- /dev/null +++ b/tests/cases/cep/spiral_TTP_2d/mesh/domain_ids.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:b7965457ae9a682c6f2599a997818ec50a7f4bbf72518b9a559f43a9b59da3c3 +size 20000 diff --git a/tests/cases/cep/spiral_TTP_2d/mesh/mesh-complete.mesh.vtu b/tests/cases/cep/spiral_TTP_2d/mesh/mesh-complete.mesh.vtu new file mode 100644 index 000000000..1f398a0e2 --- /dev/null +++ b/tests/cases/cep/spiral_TTP_2d/mesh/mesh-complete.mesh.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:e6d068e474c097975a236937fc5d3f94713456c78466c859ab726b9e95c5751a +size 164951 diff --git a/tests/cases/cep/spiral_TTP_2d/result_001.vtu b/tests/cases/cep/spiral_TTP_2d/result_001.vtu new file mode 100644 index 000000000..973f62f0b --- /dev/null +++ b/tests/cases/cep/spiral_TTP_2d/result_001.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1942cb07e1b03acafd22868527a4d44a8dcb363f97587a2cc559eb40b4dba3d5 +size 147975 diff --git a/tests/cases/cep/spiral_TTP_2d/solver.xml b/tests/cases/cep/spiral_TTP_2d/solver.xml new file mode 100644 index 000000000..37dd6ba12 --- /dev/null +++ b/tests/cases/cep/spiral_TTP_2d/solver.xml @@ -0,0 +1,81 @@ + + + + + false + 2 + 1 + 0.1 + 0.50 + STOP_SIM + + true + result + 1 + 1 + + 1 + 0 + + 1 + 0 + 0 + + + + + mesh/mesh-complete.mesh.vtu + + mesh/domain_ids.dat + + 100 + + + + + true + 1 + 3 + 1e-11 + + + TTP + domain-1_EPmodel.xml + RK + + + + TTP + domain-2_EPmodel.xml + RK + + + + TTP + domain-3_EPmodel.xml + RK + + + + TTP + domain-4_EPmodel.xml + RK + + + + true + + + + + fsils + + 500 + 1e-12 + 50 + + + + + + diff --git a/tests/cases/cep/ttp_parameters/ttp_endocardium_parameters.xml b/tests/cases/cep/ttp_parameters/ttp_endocardium_parameters.xml new file mode 100644 index 000000000..d5bdc655b --- /dev/null +++ b/tests/cases/cep/ttp_parameters/ttp_endocardium_parameters.xml @@ -0,0 +1,36 @@ + + + + 0.392 + 0.073 + + endocardium + + + + -86.709 + 138.4 + 10.355 + 1.3E-4 + 3.6E-4 + 3.715 + 0.9068 + + + + 4.48E-3 + 0.476 + 8.7E-3 + 1.55E-3 + 0.7573 + 0.7225 + 3.164E-5 + 0.8009 + 0.9778 + 0.9953 + 0.3212 + 2.235E-8 + + + + \ No newline at end of file diff --git a/tests/cases/cep/ttp_parameters/ttp_epicardium_parameters.xml b/tests/cases/cep/ttp_parameters/ttp_epicardium_parameters.xml new file mode 100644 index 000000000..d7fd60a12 --- /dev/null +++ b/tests/cases/cep/ttp_parameters/ttp_epicardium_parameters.xml @@ -0,0 +1,36 @@ + + + + 0.392 + 0.294 + + epicardium + + + + -85.23 + 136.89 + 8.6040 + 1.26E-4 + 3.6E-4 + 3.64 + 0.9073 + + + + 6.21E-3 + 0.4712 + 9.5E-3 + 1.72E-3 + 0.7444 + 0.7045 + 3.373E-5 + 0.7888 + 0.9755 + 0.9953 + 0.999998 + 2.42E-8 + + + + \ No newline at end of file diff --git a/tests/cases/cep/ttp_parameters/ttp_midmyocardium_parameters.xml b/tests/cases/cep/ttp_parameters/ttp_midmyocardium_parameters.xml new file mode 100644 index 000000000..ffd6655ff --- /dev/null +++ b/tests/cases/cep/ttp_parameters/ttp_midmyocardium_parameters.xml @@ -0,0 +1,36 @@ + + + + 0.098 + 0.294 + + myocardium + + + + -85.423 + 138.52 + 10.132 + 1.53E-4 + 4.2E-4 + 4.272 + 0.8978 + + + + 1.65E-2 + 0.4730 + 1.74E-2 + 1.65E-3 + 0.7490 + 0.6788 + 3.288E-5 + 0.7026 + 0.9526 + 0.9942 + 0.999998 + 2.347E-8 + + + + \ No newline at end of file diff --git a/tests/test_cep.py b/tests/test_cep.py index f0e49c8f1..5e6cdccba 100644 --- a/tests/test_cep.py +++ b/tests/test_cep.py @@ -22,6 +22,11 @@ def test_spiral_BO_2d(n_proc): run_with_reference(base_folder, test_folder, fields, n_proc) +def test_spiral_TTP_2d(n_proc): + test_folder = "spiral_TTP_2d" + run_with_reference(base_folder, test_folder, fields, n_proc) + + def test_square_AP_2d(n_proc): test_folder = "square_AP_2d" run_with_reference(base_folder, test_folder, fields, n_proc)