From 9259a6d28ae0e31f86dad2f231fa900ae4b205da Mon Sep 17 00:00:00 2001 From: KB Date: Thu, 23 Jul 2026 15:05:01 -0700 Subject: [PATCH 01/13] add land-niederer model --- Code/Source/solver/ActiveStress.cpp | 13 ++ Code/Source/solver/ActiveStress.h | 11 + Code/Source/solver/CMakeLists.txt | 1 + Code/Source/solver/Parameters.cpp | 9 + Code/Source/solver/Parameters.h | 6 + .../solver/active_stress_land_niederer.cpp | 157 ++++++++++++++ .../solver/active_stress_land_niederer.h | 143 +++++++++++++ .../slab_LandNiederer/mesh/X0.vtp | 3 + .../slab_LandNiederer/mesh/X1.vtp | 3 + .../slab_LandNiederer/mesh/volume.vtu | 3 + .../slab_LandNiederer/solver.xml | 199 ++++++++++++++++++ 11 files changed, 548 insertions(+) create mode 100644 Code/Source/solver/active_stress_land_niederer.cpp create mode 100644 Code/Source/solver/active_stress_land_niederer.h create mode 100644 tests/cases/electromechanics/slab_LandNiederer/mesh/X0.vtp create mode 100644 tests/cases/electromechanics/slab_LandNiederer/mesh/X1.vtp create mode 100644 tests/cases/electromechanics/slab_LandNiederer/mesh/volume.vtu create mode 100644 tests/cases/electromechanics/slab_LandNiederer/solver.xml diff --git a/Code/Source/solver/ActiveStress.cpp b/Code/Source/solver/ActiveStress.cpp index 7db495dae..32084889f 100644 --- a/Code/Source/solver/ActiveStress.cpp +++ b/Code/Source/solver/ActiveStress.cpp @@ -14,6 +14,8 @@ void ActiveStress::read_parameters(const ActiveStressParameters ¶ms) { eta_s = params.get_eta_s(); eta_n = params.get_eta_n(); + use_stabilization = params.get_use_stabilization(); + read_model_specific_parameters( params.get_parameters(params.get_model_name())); } @@ -23,6 +25,7 @@ void ActiveStress::distribute_parameters(const CmMod &cm_mod, cm.bcast(cm_mod, &eta_f); cm.bcast(cm_mod, &eta_s); cm.bcast(cm_mod, &eta_n); + cm.bcast(cm_mod, &use_stabilization); distribute_model_specific_parameters(cm_mod, cm); } @@ -40,6 +43,16 @@ void ActiveStress::init(const unsigned int tnNo) { } active_tension.resize(tnNo); + raw_active_tension.resize(tnNo); + previous_fiber_stretch.resize(tnNo); + has_previous_fiber_stretch.resize(tnNo); + + for (unsigned int i = 0; i < tnNo; ++i) { + active_tension[i] = 0.0; + raw_active_tension[i] = 0.0; + previous_fiber_stretch[i] = 1.0; + has_previous_fiber_stretch[i] = 0; + } } void ActiveStress::advance_time_step(const double t, const double dt, diff --git a/Code/Source/solver/ActiveStress.h b/Code/Source/solver/ActiveStress.h index 29a419cf6..fea5f31c4 100644 --- a/Code/Source/solver/ActiveStress.h +++ b/Code/Source/solver/ActiveStress.h @@ -275,6 +275,17 @@ class ActiveStress { /// Active tension at every node. Vector active_tension; + /// Raw active tension at every node, before stabilization is applied. + Vector raw_active_tension; + + /// Previous fiber stretch at every node, used for stabilization. + Vector previous_fiber_stretch; + + /// Whether to apply stabilization to the active tension. + bool use_stabilization = false; + + Vector has_previous_fiber_stretch; + /// Active tension coefficient along the fiber direction. double eta_f; diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt index 468edf0fd..7f8d1d775 100644 --- a/Code/Source/solver/CMakeLists.txt +++ b/Code/Source/solver/CMakeLists.txt @@ -269,6 +269,7 @@ set(CSRCS ActiveStressODE.cpp ActiveStressNashPanfilov.cpp ActiveStressRegazzoni.cpp + active_stress_land_niederer.cpp SPLIT.c diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index f0922a258..5d84c260e 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -1910,7 +1910,12 @@ const std::string ActiveStressParameters::xml_element_name = "Active_stress"; ActiveStressParameters::ActiveStressParameters() { model_name = Parameter("Model", "", true); + use_stabilization = + Parameter("Use_stabilization", false, false); + set_parameter("Model", "", /* required = */ true, model_name); + set_parameter("Use_stabilization", false, /* required = */ false, + use_stabilization); ActiveStressFactory::visit( [this](const std::string &name, const ActiveStress &model) { @@ -1984,6 +1989,10 @@ double ActiveStressParameters::get_eta_n() const { return directional_distribution.sheet_normal_direction.value(); } +bool ActiveStressParameters::get_use_stabilization() const { + return use_stabilization.value(); +} + const ActiveStressModelParameters & ActiveStressParameters::get_parameters(const std::string &model_name) const { if (active_stress_models.count(model_name) == 0) { diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index 94d3b677b..aa3e7a627 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1534,6 +1534,9 @@ class ActiveStressParameters : public ParameterLists { /// Get the active tension coefficient along sheet normals. double get_eta_n() const; + /// Get the stabilization flag for Regazzoni-style active tension stabilization. + bool get_use_stabilization() const; + /// Get the parameters for a given active stress model. const ActiveStressModelParameters & get_parameters(const std::string &model_name) const; @@ -1545,6 +1548,9 @@ class ActiveStressParameters : public ParameterLists { /// Parameter for the model name. Parameter model_name; + /// Flag for Regazzoni-style active tension stabilization + Parameter use_stabilization; + /// Parameters for the directional distribution of active tension. DirectionalDistributionParameters directional_distribution; diff --git a/Code/Source/solver/active_stress_land_niederer.cpp b/Code/Source/solver/active_stress_land_niederer.cpp new file mode 100644 index 000000000..440476e7e --- /dev/null +++ b/Code/Source/solver/active_stress_land_niederer.cpp @@ -0,0 +1,157 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause + +#include "active_stress_land_niederer.h" + +void LandNiederer::read_model_specific_parameters( + const ActiveStressModelParameters ¶ms) { + ActiveStressODE::read_model_specific_parameters(params); + + CaRef = params.get_scalar("CaRef"); + eta_Tm = params.get_scalar("eta_Tm"); + k_uw = params.get_scalar("k_uw"); + k_ws = params.get_scalar("k_ws"); + Tref = params.get_scalar("Tref"); + k_TRPN = params.get_scalar("k_TRPN"); + eta_TRPN = params.get_scalar("eta_TRPN"); + k_u = params.get_scalar("k_u"); + TRPN50 = params.get_scalar("TRPN50"); + rw = params.get_scalar("rw"); + rs = params.get_scalar("rs"); + gamma_s = params.get_scalar("gamma_s"); + gamma_w = params.get_scalar("gamma_w"); + phi = params.get_scalar("phi"); + Aeff = params.get_scalar("Aeff"); + beta0 = params.get_scalar("beta0"); + beta1 = params.get_scalar("beta1"); +} + +void LandNiederer::distribute_model_specific_parameters(const CmMod &cm_mod, + const cmType &cm) { + ActiveStressODE::distribute_model_specific_parameters(cm_mod, cm); + + cm.bcast(cm_mod, &CaRef); + cm.bcast(cm_mod, &eta_Tm); + cm.bcast(cm_mod, &k_uw); + cm.bcast(cm_mod, &k_ws); + cm.bcast(cm_mod, &Tref); + cm.bcast(cm_mod, &k_TRPN); + cm.bcast(cm_mod, &eta_TRPN); + cm.bcast(cm_mod, &k_u); + cm.bcast(cm_mod, &TRPN50); + cm.bcast(cm_mod, &rw); + cm.bcast(cm_mod, &rs); + cm.bcast(cm_mod, &gamma_s); + cm.bcast(cm_mod, &gamma_w); + cm.bcast(cm_mod, &phi); + cm.bcast(cm_mod, &Aeff); + cm.bcast(cm_mod, &beta0); + cm.bcast(cm_mod, &beta1); +} + +void LandNiederer::init_local(Vector &state) const { + state[0] = 1.0; + state[1] = 0.0; + state[2] = 0.0; + state[3] = 0.0; + state[4] = 0.0; + state[5] = 0.0; + const double lambda = 1.0; + state[6] = std::max(0.0, 1.0 + beta0 * (lambda + std::min(0.87, lambda) - 1.87)); +} + +void LandNiederer::advance_time_step_local(const double t, const double dt, + const double calcium, + const double fiber_stretch, + const double fiber_stretch_rate, + Vector &state) const { + + Vector ode_state(6); + for (int i = 0; i < 6; i++) { + ode_state[i] = state[i]; + } + + ActiveStressODE::advance_time_step_local(t, dt, calcium, fiber_stretch, + fiber_stretch_rate, ode_state); + + for (unsigned int i = 0; i < 6; ++i) { + state[i] = ode_state[i]; + } + + const double lambda = std::min(1.2, fiber_stretch); + state[6] = std::max(0.0, 1.0 + beta0 * (lambda + std::min(0.87, lambda) - 1.87)); +} + +Vector LandNiederer::getf(const double t, const Vector &state, + const double calcium, + const double fiber_stretch, + const double fiber_stretch_rate) const { + Vector f(6); + + // State Variables + const double XB = state[0]; + const double XW = std::max(0.0, state[1]); + const double CaTRPN = std::max(0.0, state[2]); + const double XS = std::max(0.0, state[3]); + const double ZETAW = state[4]; + const double ZETAS = state[5]; + + // Bound State Equation + const double lambda = std::min(1.2, fiber_stretch); + const double XU = (1 - XB) - XS - XW; + const double k_b = k_u * std::pow(TRPN50, eta_Tm)/ (1 - rs - (1 - rs)*rw); + f[0] = k_b * std::min(100.0, std::pow(CaTRPN, -(eta_Tm/2.0))) * XU - k_u * std::pow(CaTRPN, (eta_Tm/2.0)) * XB; + + // Pre-power Stroke Kinetics + const double k_wu = (k_uw * (1.0/rw - 1.0) - k_ws); + const double gamma_wu = gamma_w * std::abs(ZETAW); + f[1] = k_uw*XU - k_wu*XW - k_ws*XW - gamma_wu*XW; + + // Calcium-Troponin Binding Kinetics + const double CaT50 = CaRef + beta1*std::min(0.2, lambda - 1.0); + f[2] = k_TRPN * (std::pow((calcium/CaT50), eta_TRPN) * (1.0 - CaTRPN) - CaTRPN); + + // Post-power Stroke Kinetics + const double k_su = k_ws * rw * (1.0/rs - 1.0); + const double term1 = (ZETAS > 0.0) ? ZETAS : 0.0; + const double term2 = (ZETAS < -1.0) ? (-ZETAS - 1.0) : 0.0; + const double gamma_su = gamma_s * std::max(term1, term2); + f[3] = k_ws*XW - k_su*XS - gamma_su*XS; + + // Distortion-Decay Kinetics + const double Aw = Aeff * rs/((1.0 - rs) * rw + rs); + const double cw = phi * k_uw * (1.0 - rs)*(1.0 - rw) / ((1.0 - rs)*rw); + f[4] = Aw * fiber_stretch_rate - cw*ZETAW; + + const double As = Aw; + const double cs = phi * k_ws * (1.0 - rs)*rw/rs; + f[5] = As * fiber_stretch_rate - cs*ZETAS; + + return f; +} + +double LandNiederer::compute_active_tension_local(const Vector &state) const { + const double XW = std::max(0.0, state[1]); + const double XS = std::max(0.0, state[3]); + const double ZETAW = state[4]; + const double ZETAS = state[5]; + const double LFac = state[6]; + + return LFac * (Tref/rs) * ((ZETAS+1.0) * XS + (ZETAW) * XW); +} + + +double LandNiederer::compute_active_stiffness_local(const Vector &state, + const double fiber_stretch, + const double fiber_stretch_rate) const { + const double XW = std::max(0.0, state[1]); + const double XS = std::max(0.0, state[3]); + const double LFac = state[6]; + const double Aw = Aeff * rs/((1.0 - rs) * rw + rs); + const double As = Aw; + + return LFac * (Tref/rs) * (As*XS + Aw*XW); +} + + +REGISTER_ACTIVE_STRESS_MODEL("LandNiederer", LandNiederer); \ No newline at end of file diff --git a/Code/Source/solver/active_stress_land_niederer.h b/Code/Source/solver/active_stress_land_niederer.h new file mode 100644 index 000000000..4d08c4ac2 --- /dev/null +++ b/Code/Source/solver/active_stress_land_niederer.h @@ -0,0 +1,143 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause + +#ifndef ACTIVE_STRESS_LAND_NIEDERER_H +#define ACTIVE_STRESS_LAND_NIEDERER_H + +#include "ActiveStressODE.h" + +class LandNiederer : public ActiveStressODE { +public: + /// Model label. + static inline const std::string label = "LandNiederer"; + + /// Model parameters class. + class Parameters : public ActiveStressODE::Parameters { + public: + Parameters() : ActiveStressODE::Parameters(label) { + constexpr bool required = true; + + add_parameter("CaRef", 0.805, required); + add_parameter("eta_Tm", 5.0, required); + add_parameter("k_uw", 0.182, required); + add_parameter("k_ws", 0.012, required); + add_parameter("Tref", 120.0, required); + add_parameter("k_TRPN", 0.1, required); + add_parameter("eta_TRPN", 2.0, required); + add_parameter("k_u", 1.0, required); + add_parameter("TRPN50", 0.35, required); + add_parameter("rw", 0.5, required); + add_parameter("rs", 0.25, required); + add_parameter("gamma_s", 0.0085, required); + add_parameter("gamma_w", 0.615, required); + add_parameter("phi", 2.23, required); + add_parameter("Aeff", 25.0, required); + add_parameter("beta0", 2.3, required); + add_parameter("beta1", -2.4, required); + } + }; + + /** + * @brief Constructor. + */ + LandNiederer() : ActiveStressODE(7) {} + + /** + * @brief Construct an instance of model parameters. + */ + virtual std::unique_ptr + get_parameters() const override { + return std::make_unique(); + } + +protected: + /** + * @brief Read model parameters from a parameter object. + */ + virtual void read_model_specific_parameters( + const ActiveStressModelParameters ¶ms) override; + + /** + * @brief Distribute model parameters to all parallel processes. + */ + virtual void distribute_model_specific_parameters(const CmMod &cm_mod, + const cmType &cm) override; + + /** + * @brief Initialize the state vector for a single node. + * + * @param[out] state State vector for a single node, to be initialized by + * this function. + */ + virtual void init_local(Vector &state) const override; + + virtual void advance_time_step_local(const double t, const double dt, + const double calcium, + const double fiber_stretch, + const double fiber_stretch_rate, + Vector &state) const override; + + /** + * @brief Compute the rate of change in the state variables. + */ + virtual Vector getf(const double t, const Vector &state, + const double calcium, const double fiber_stretch, + const double fiber_stretch_rate) const override; + + /** + * @brief Compute the active tension for a single node. + */ + virtual double + compute_active_tension_local(const Vector &state) const override; + + /** + * @brief Compute the active stiffness for a single node. + */ + virtual double + compute_active_stiffness_local(const Vector &state, + const double fiber_stretch, + const double fiber_stretch_rate) const override; + + /// @name Model parameters. + /// @{ + + double CaRef; + + double eta_Tm; + + double k_uw; + + double k_ws; + + double Tref; + + double k_TRPN; + + double eta_TRPN; + + double k_u; + + double TRPN50; + + double rw; + + double rs; + + double gamma_s; + + double gamma_w; + + double phi; + + double Aeff; + + double beta0; + + double beta1; + + /// @} +}; + + + +#endif \ No newline at end of file diff --git a/tests/cases/electromechanics/slab_LandNiederer/mesh/X0.vtp b/tests/cases/electromechanics/slab_LandNiederer/mesh/X0.vtp new file mode 100644 index 000000000..eaebcc0eb --- /dev/null +++ b/tests/cases/electromechanics/slab_LandNiederer/mesh/X0.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:3c0889ad4a8a4659996309b5bba825c1ac5c3f03d2faca2e84e48f1e035ca240 +size 5219 diff --git a/tests/cases/electromechanics/slab_LandNiederer/mesh/X1.vtp b/tests/cases/electromechanics/slab_LandNiederer/mesh/X1.vtp new file mode 100644 index 000000000..07ad82435 --- /dev/null +++ b/tests/cases/electromechanics/slab_LandNiederer/mesh/X1.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:8248936f54f9537dcad8e2344beb804485d0f5cbfb960f0a8b025310e1af99ae +size 5128 diff --git a/tests/cases/electromechanics/slab_LandNiederer/mesh/volume.vtu b/tests/cases/electromechanics/slab_LandNiederer/mesh/volume.vtu new file mode 100644 index 000000000..dd12aca89 --- /dev/null +++ b/tests/cases/electromechanics/slab_LandNiederer/mesh/volume.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:2d9c504f6d9dc221739d107e19414233cdda640d73647eba03902a1605274149 +size 421967 diff --git a/tests/cases/electromechanics/slab_LandNiederer/solver.xml b/tests/cases/electromechanics/slab_LandNiederer/solver.xml new file mode 100644 index 000000000..e8f08e9da --- /dev/null +++ b/tests/cases/electromechanics/slab_LandNiederer/solver.xml @@ -0,0 +1,199 @@ + + + + false + 3 + 100 + 1.0 + 0.50 + STOP_SIM + + true + result + 10 + 0 + + 1000 + 0 + + 1 + 1 + 0 + + + + ./mesh/volume.vtu + + + ./mesh/X0.vtp + + + + ./mesh/X1.vtp + + + ./mesh/volume.vtu + + (1, 0, 0) + (0, 1, 0) + (0, 1, 0) + + + + true + + 1 + 1 + 1e-6 + + + TTP + + 0.012571 + 0.082715 + 0.0 + 0.0 + + ../../cep/ttp_parameters/ttp_epicardium_parameters.xml + + + 14.838 + 3.98E-5 + 0.153 + + + RK4 + + + + TTP + + 0.012571 + 0.082715 + 0.0 + 0.0 + + ../../cep/ttp_parameters/ttp_epicardium_parameters.xml + + + 14.838 + 3.98E-5 + 0.153 + + + + -35.714 + 0.0 + 2.0 + 10000.0 + + + RK4 + + + + true + true + + + + + fsils + + 100 + 1e-6 + 50 + + + + + 1 + 6 + 1e-3 + + 1e-3 + + + 59.0e-6 + 8.023 + 18472.0e-6 + 16.026 + 2481.0e-6 + 11.12 + 216.0e-6 + 11.436 + 100.0 + + + ST91 + 1.0 + + + 1.0 + + + + LandNiederer + true + + + 0.7 + 0.2 + 0.1 + + + + FE + + 0.805 + 5.0 + 0.182 + 0.012 + 120.0 + 0.1 + 2.0 + 1.0 + 0.35 + 0.5 + 0.25 + 0.0085 + 0.615 + 2.23 + 25.0 + 2.3 + -2.4 + + + + + true + true + true + true + true + true + true + true + + true + true + true + + + + + fsils + + 1e-4 + 1e-6 + 1000 + + + + Dir + 0.0 + + + + + From 0e8f0cff4db9ac2173aec5de95e0f12deca4ac4f Mon Sep 17 00:00:00 2001 From: KB Date: Mon, 27 Jul 2026 15:56:55 -0700 Subject: [PATCH 02/13] unit test and test case update --- .../solver/active_stress_land_niederer.cpp | 5 +- .../solver/active_stress_land_niederer.h | 20 +++- .../slab_LandNiederer/solver.xml | 48 ++++---- .../Testing/Temporary/CTestCostData.txt | 1 + .../unitTests/Testing/Temporary/LastTest.log | 3 + .../test_electromechanics_models.cpp | 104 ++++++++++++++++++ 6 files changed, 153 insertions(+), 28 deletions(-) create mode 100644 tests/unitTests/Testing/Temporary/CTestCostData.txt create mode 100644 tests/unitTests/Testing/Temporary/LastTest.log create mode 100644 tests/unitTests/ionic_model_tests/test_electromechanics_models.cpp diff --git a/Code/Source/solver/active_stress_land_niederer.cpp b/Code/Source/solver/active_stress_land_niederer.cpp index 440476e7e..4590ee18c 100644 --- a/Code/Source/solver/active_stress_land_niederer.cpp +++ b/Code/Source/solver/active_stress_land_niederer.cpp @@ -71,8 +71,11 @@ void LandNiederer::advance_time_step_local(const double t, const double dt, ode_state[i] = state[i]; } + /// TODO: Remove the forced fiber_stretch_rate = 0.0 below and use the actual fiber_stretch_rate. + /// This is a temporary fix to avoid the active stress model from blowing up when the fiber_stretch_rate is large. + ActiveStressODE::advance_time_step_local(t, dt, calcium, fiber_stretch, - fiber_stretch_rate, ode_state); + 0.0, ode_state); for (unsigned int i = 0; i < 6; ++i) { state[i] = ode_state[i]; diff --git a/Code/Source/solver/active_stress_land_niederer.h b/Code/Source/solver/active_stress_land_niederer.h index 4d08c4ac2..9c44ae861 100644 --- a/Code/Source/solver/active_stress_land_niederer.h +++ b/Code/Source/solver/active_stress_land_niederer.h @@ -100,42 +100,56 @@ class LandNiederer : public ActiveStressODE { /// @name Model parameters. /// @{ - + /// Calcium Sensitivity [uM] double CaRef; + /// Cooperativity coefficient [-] double eta_Tm; + /// Rate constant for transition from unbound (U) to pre-power stroke (W) state [1/ms] double k_uw; + /// Rate constant for transition from pre-power stroke (W) to post-power stroke (S) state [1/ms] double k_ws; + /// Maximum observable tension @f$\Tref@f$ at resting length [kPa] double Tref; + /// Unbinding rate constant for calcium-troponin C complex [1/ms] double k_TRPN; + /// Cooperativity of the calcium-troponin C binding rate [-] double eta_TRPN; + /// Unbinding rate constant [1/ms] double k_u; + /// The value of CaTRPN where bounded state B = 0.5 in steady-state [-] double TRPN50; + /// Steady-state ratio between pre-powerstroke and non-strongly bound [-] double rw; + /// Steady-state duty ratio [-] double rs; + /// [-] double gamma_s; + /// [-] double gamma_w; + /// [-] double phi; + /// Rescaled double Aeff; + /// Change in maximal tension based on changes in filament overlap double beta0; + /// Change in calcium sensitivity double beta1; - - /// @} }; diff --git a/tests/cases/electromechanics/slab_LandNiederer/solver.xml b/tests/cases/electromechanics/slab_LandNiederer/solver.xml index e8f08e9da..ccf80b824 100644 --- a/tests/cases/electromechanics/slab_LandNiederer/solver.xml +++ b/tests/cases/electromechanics/slab_LandNiederer/solver.xml @@ -3,7 +3,7 @@ false 3 - 100 + 10 1.0 0.50 STOP_SIM @@ -44,7 +44,7 @@ 1 1 - 1e-6 + 1e-12 TTP @@ -101,7 +101,7 @@ fsils 100 - 1e-6 + 1e-12 50 @@ -109,7 +109,7 @@ 1 6 - 1e-3 + 1e-12 1e-3 @@ -134,7 +134,7 @@ LandNiederer - true + false 0.7 @@ -145,23 +145,23 @@ FE - 0.805 - 5.0 - 0.182 - 0.012 - 120.0 - 0.1 - 2.0 - 1.0 - 0.35 - 0.5 - 0.25 - 0.0085 - 0.615 - 2.23 - 25.0 - 2.3 - -2.4 + 0.000805 + 5.0 + 0.182 + 0.012 + 0.12 + 0.1 + 2.0 + 1.0 + 0.35 + 0.5 + 0.25 + 0.0085 + 0.615 + 2.23 + 25.0 + 0.0 + 0.0 @@ -184,8 +184,8 @@ fsils - 1e-4 - 1e-6 + 1e-12 + 1e-14 1000 diff --git a/tests/unitTests/Testing/Temporary/CTestCostData.txt b/tests/unitTests/Testing/Temporary/CTestCostData.txt new file mode 100644 index 000000000..ed97d539c --- /dev/null +++ b/tests/unitTests/Testing/Temporary/CTestCostData.txt @@ -0,0 +1 @@ +--- diff --git a/tests/unitTests/Testing/Temporary/LastTest.log b/tests/unitTests/Testing/Temporary/LastTest.log new file mode 100644 index 000000000..d25f1e60c --- /dev/null +++ b/tests/unitTests/Testing/Temporary/LastTest.log @@ -0,0 +1,3 @@ +Start testing: Jul 27 11:37 PDT +---------------------------------------------------------- +End testing: Jul 27 11:37 PDT diff --git a/tests/unitTests/ionic_model_tests/test_electromechanics_models.cpp b/tests/unitTests/ionic_model_tests/test_electromechanics_models.cpp new file mode 100644 index 000000000..941e34a70 --- /dev/null +++ b/tests/unitTests/ionic_model_tests/test_electromechanics_models.cpp @@ -0,0 +1,104 @@ +/* Copyright (c) Stanford University, The Regents of the University of + * California, and others. + * + * All Rights Reserved. + * + * See Copyright-SimVascular.txt for additional details. + * + * Permission is hereby granted, free of charge, to any person obtaining + * a copy of this software and associated documentation files (the + * "Software"), to deal in the Software without restriction, including + * without limitation the rights to use, copy, modify, merge, publish, + * distribute, sublicense, and/or sell copies of the Software, and to + * permit persons to whom the Software is furnished to do so, subject + * to the following conditions: + * + * The above copyright notice and this permission notice shall be included + * in all copies or substantial portions of the Software. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS + * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED + * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A + * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER + * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, + * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, + * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR + * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF + * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING + * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + */ + +#include "../test_common.h" +#include "active_stress_land_niederer.h" + +class TestableLandNiederer : public LandNiederer { +public: + Vector evaluate_getf(const double t, const Vector &state, + const double calcium, + const double fiber_stretch, + const double fiber_stretch_rate) const { + return getf(t, state, calcium, fiber_stretch, fiber_stretch_rate); + } + + void set_parameters_for_test() { + CaRef = 0.805; + eta_Tm = 5.0; + k_uw = 0.182; + k_ws = 0.012; + Tref = 120.0; + k_TRPN = 0.1; + eta_TRPN = 2.0; + k_u = 1.0; + TRPN50 = 0.35; + rw = 0.5; + rs = 0.25; + gamma_s = 0.0085; + gamma_w = 0.615; + phi = 2.23; + Aeff = 25.0; + beta0 = 2.3; + beta1 = -2.4; + } +}; + +class ElectromechanicsModelTest : public ::testing::Test { +protected: + void SetUp() override { model.set_parameters_for_test(); } + + void TearDown() override {} + + static Vector initial_state() { + return {0.8900, 0.0444, 0.2609, 0.0052, 0.0000, 0.0000}; + } + + static Vector expected_rhs() { + return {-0.0066, 0.0029, 0.0029, 0.0004, 1.0, 1.0}; + } + + TestableLandNiederer model; +}; + +TEST_F(ElectromechanicsModelTest, LandNiedererGetf) { + // Evaluate the Land-Niederer ODE right-hand side at one state and compare + // against reference values computed from the MATLAB implementation. + + const double t = 0.0; + const double calcium = 0.5040; + const double fiber_stretch = 1.0; + const double fiber_stretch_rate = 0.1; + + const Vector state = initial_state(); + const Vector expected = expected_rhs(); + + const Vector rhs = + model.evaluate_getf(t, state, calcium, fiber_stretch, + fiber_stretch_rate); + + ASSERT_EQ(rhs.size(), expected.size()); + + for (int i = 0; i < rhs.size(); ++i) { + ASSERT_NEAR(rhs[i], expected[i], 1e-3) + << "RHS mismatch at component " << i; + } +} \ No newline at end of file From 3176722151b2fd100c2e6493dfa7e35ed577f063 Mon Sep 17 00:00:00 2001 From: KB Date: Mon, 3 Aug 2026 08:00:01 -0700 Subject: [PATCH 03/13] Added RK4 to stabilize Land Niederer ODE --- Code/Source/solver/ActiveStressODE.cpp | 15 ++++++++ Code/Source/solver/ActiveStressODE.h | 37 ++++++++++++++++--- .../solver/active_stress_land_niederer.cpp | 5 +-- 3 files changed, 47 insertions(+), 10 deletions(-) diff --git a/Code/Source/solver/ActiveStressODE.cpp b/Code/Source/solver/ActiveStressODE.cpp index 9b9ca9e17..6fa2b92b1 100644 --- a/Code/Source/solver/ActiveStressODE.cpp +++ b/Code/Source/solver/ActiveStressODE.cpp @@ -9,6 +9,8 @@ void ActiveStressODE::read_model_specific_parameters( if (solver_str == "FE") { ode_solver = ODESolver::ForwardEuler; + } else if (solver_str == "RK4" || solver_str == "RK") { + ode_solver = ODESolver::RungeKutta4; } else { svmp::raise("Unknown ODE solver " + solver_str + " for active stress models."); @@ -29,6 +31,19 @@ void ActiveStressODE::advance_time_step_local(const double t, const double dt, const Vector f = getf(t - dt, state, calcium, fiber_stretch, fiber_stretch_rate); state.add(dt, f); + } else if (ode_solver == ODESolver::RungeKutta4) { + const Vector k1 = + getf(t - dt, state, calcium, fiber_stretch, fiber_stretch_rate); + const Vector k2 = + getf(t - 0.5 * dt, state + 0.5 * dt * k1, calcium, fiber_stretch, + fiber_stretch_rate); + const Vector k3 = + getf(t - 0.5 * dt, state + 0.5 * dt * k2, calcium, fiber_stretch, + fiber_stretch_rate); + const Vector k4 = getf(t, state + dt * k3, calcium, fiber_stretch, + fiber_stretch_rate); + + state.add(dt / 6.0, k1 + 2.0 * (k2 + k3) + k4); } else { svmp::raise( "Unknown ODE solver " + std::to_string(static_cast(ode_solver)) + diff --git a/Code/Source/solver/ActiveStressODE.h b/Code/Source/solver/ActiveStressODE.h index 3efd47a31..6b573d71e 100644 --- a/Code/Source/solver/ActiveStressODE.h +++ b/Code/Source/solver/ActiveStressODE.h @@ -69,7 +69,34 @@ class ActiveStressODE : public ActiveStress { * \fiberstretchrate_i^n)\;. * @f] */ - ForwardEuler + ForwardEuler, + + /** + * @brief 4th order explicit Runge-Kutta. + * + * The state vector is updated as follows: + * @f[ \begin{aligned} + * \mathbf{k}_1 &= \mathbf{F}_\text{AS}(t^n, \astressstate_i^n, + * \calcium_i^n, \fiberstretch_i^n, \fiberstretchrate_i^n)\;, \\ + * \mathbf{k}_2 &= \mathbf{F}_\text{AS}(t^n + \Delta t / 2, + * \astressstate_i^n + \Delta t \mathbf{k}_1 / 2, \calcium_i^n, + * \fiberstretch_i^n, \fiberstretchrate_i^n)\;, \\ + * \mathbf{k}_3 &= \mathbf{F}_\text{AS}(t^n + \Delta t / 2, + * \astressstate_i^n + \Delta t \mathbf{k}_2 / 2, \calcium_i^n, + * \fiberstretch_i^n, \fiberstretchrate_i^n)\;, \\ + * \mathbf{k}_4 &= \mathbf{F}_\text{AS}(t^{n+1}, + * \astressstate_i^n + \Delta t \mathbf{k}_3, \calcium_i^n, + * \fiberstretch_i^n, \fiberstretchrate_i^n)\;, \\ + * \astressstate_i^{n+1} &= \astressstate_i^n + \frac{\Delta t}{6} + * \left( \mathbf{k}_1 + 2 \mathbf{k}_2 + 2 \mathbf{k}_3 + * + \mathbf{k}_4 \right)\;. + * \end{aligned} @f] + * + * Note that the driving quantities @f$\calcium@f$, @f$\fiberstretch@f$ and + * @f$\fiberstretchrate@f$ are held fixed at their values at @f$t^{n+1}@f$ + * across all four stages. + */ + RungeKutta4 }; /** @@ -125,11 +152,9 @@ class ActiveStressODE : public ActiveStress { * @param[in,out] state State vector for a single node, to be updated by * this function. * - * @todo[michelebucelli] It might be necessary or useful to implement other - * timestepping schemes, e.g. Runge-Kutta. In that case, we might want to - * expand the interface to support implicit time stepping too, e.g. by - * adding a method to evaluate the Jacobian matrix of the system, as in - * @ref IonicModel. + * @todo[michelebucelli] It might be necessary or useful to implement + * implicit time stepping too, e.g. by adding a method to evaluate the + * Jacobian matrix of the system, as in @ref IonicModel. */ virtual void advance_time_step_local(const double t, const double dt, const double calcium, diff --git a/Code/Source/solver/active_stress_land_niederer.cpp b/Code/Source/solver/active_stress_land_niederer.cpp index 4590ee18c..440476e7e 100644 --- a/Code/Source/solver/active_stress_land_niederer.cpp +++ b/Code/Source/solver/active_stress_land_niederer.cpp @@ -71,11 +71,8 @@ void LandNiederer::advance_time_step_local(const double t, const double dt, ode_state[i] = state[i]; } - /// TODO: Remove the forced fiber_stretch_rate = 0.0 below and use the actual fiber_stretch_rate. - /// This is a temporary fix to avoid the active stress model from blowing up when the fiber_stretch_rate is large. - ActiveStressODE::advance_time_step_local(t, dt, calcium, fiber_stretch, - 0.0, ode_state); + fiber_stretch_rate, ode_state); for (unsigned int i = 0; i < 6; ++i) { state[i] = ode_state[i]; From ed396a9c18b6deec0d1783186ac75a444baa24b5 Mon Sep 17 00:00:00 2001 From: KB Date: Wed, 19 Aug 2026 22:07:59 -0700 Subject: [PATCH 04/13] added more comments - will remove some verbose parts later --- .../solver/active_stress_land_niederer.h | 103 +++++++++++------- 1 file changed, 65 insertions(+), 38 deletions(-) diff --git a/Code/Source/solver/active_stress_land_niederer.h b/Code/Source/solver/active_stress_land_niederer.h index 9c44ae861..669bca775 100644 --- a/Code/Source/solver/active_stress_land_niederer.h +++ b/Code/Source/solver/active_stress_land_niederer.h @@ -8,15 +8,20 @@ class LandNiederer : public ActiveStressODE { public: - /// Model label. + /// Model label, used for factory registration and XML selection. static inline const std::string label = "LandNiederer"; - /// Model parameters class. + /* + * @brief Model parameters class. + * Declares the parameters required by the model. All parameters are + * marked as required, and omitting a parameter will cause a parse error. + */ class Parameters : public ActiveStressODE::Parameters { public: Parameters() : ActiveStressODE::Parameters(label) { constexpr bool required = true; + // Reference values: Land et al. 2017 add_parameter("CaRef", 0.805, required); add_parameter("eta_Tm", 5.0, required); add_parameter("k_uw", 0.182, required); @@ -100,13 +105,17 @@ class LandNiederer : public ActiveStressODE { /// @name Model parameters. /// @{ - /// Calcium Sensitivity [uM] + /// Reference intracellular calcium concentration giving half-maximal + /// troponin C saturation, @f$[Ca^{2+}]_{T50,ref}@f$, used (with length + /// dependence via beta1) in the CaTRPN binding ODE [uM] double CaRef; - /// Cooperativity coefficient [-] + /// Cooperativity (Hill) exponent @f$n_{Tm}@f$ for the tropomyosin + /// blocked/unblocked (B/U) transition; sets the steepness of thin-filament + /// activation by CaTRPN [-] double eta_Tm; - /// Rate constant for transition from unbound (U) to pre-power stroke (W) state [1/ms] + /// Rate constant for transition from unbound (U) to pre-power stroke (W) state [1/ms] double k_uw; /// Rate constant for transition from pre-power stroke (W) to post-power stroke (S) state [1/ms] @@ -115,41 +124,59 @@ class LandNiederer : public ActiveStressODE { /// Maximum observable tension @f$\Tref@f$ at resting length [kPa] double Tref; - /// Unbinding rate constant for calcium-troponin C complex [1/ms] + /// Rate constant @f$k_{TRPN}@f$ governing calcium binding/unbinding + /// kinetics of troponin C (the "pace" of CaTRPN relaxation to its + /// steady-state value) [1/ms] double k_TRPN; - /// Cooperativity of the calcium-troponin C binding rate [-] - double eta_TRPN; - - /// Unbinding rate constant [1/ms] - double k_u; - - /// The value of CaTRPN where bounded state B = 0.5 in steady-state [-] - double TRPN50; - - /// Steady-state ratio between pre-powerstroke and non-strongly bound [-] - double rw; - - /// Steady-state duty ratio [-] - double rs; - - /// [-] - double gamma_s; - - /// [-] - double gamma_w; - - /// [-] - double phi; - - /// Rescaled - double Aeff; - - /// Change in maximal tension based on changes in filament overlap - double beta0; - - /// Change in calcium sensitivity - double beta1; + /// Cooperativity (Hill) exponent @f$n_{TRPN}@f$ of the calcium-troponin C + /// binding rate [-] + double eta_TRPN; + + /// Tropomyosin unblocking rate constant @f$k_u@f$ (rate at which blocked + /// binding sites become unblocked); together with rw, rs, and TRPN50 it + /// fixes the blocking rate kb [1/ms] + double k_u; + + /// The value of CaTRPN where bounded state B = 0.5 in steady-state [-] + double TRPN50; + + /// Steady-state ratio @f$r_w@f$ between the pre-powerstroke (W) population + /// and the non-strongly-bound (U+W) population at equilibrium; sets the + /// reverse rate constant kwu [-] + double rw; + + /// Steady-state duty ratio @f$r_s@f$: the fraction of cross-bridges in the + /// strongly-bound (post-powerstroke, S) state at equilibrium; sets the + /// reverse rate constant ksu [-] + double rs; + + /// Strain-dependent detachment-rate coefficient @f$\gamma_s@f$ for the + /// strongly-bound (post-powerstroke) S state; scales how fast cross-bridges + /// are pulled off by distortion in that state [-] + double gamma_s; + + /// Strain-dependent detachment-rate coefficient @f$\gamma_w@f$ for the + /// weakly-bound (pre-powerstroke) W state; scales how fast cross-bridges + /// are pulled off by distortion in that state [-] + double gamma_w; + + /// Proportionality factor @f$\phi@f$ relating the cross-bridge distortion + /// decay rates (cw, cs) to the forward cycling rates k_uw, k_ws [-] + double phi; + + /// Rescaled magnitude @f$A_{eff}@f$ of the immediate cross-bridge + /// distortion response to fiber stretch rate; sets the distortion + /// sensitivities As, Aw [-] + double Aeff; + + /// Coefficient @f$\beta_0@f$ controlling the length-dependence of maximal + /// tension through changes in filament overlap (Frank-Starling effect) [-] + double beta0; + + /// Coefficient @f$\beta_1@f$ controlling the length-dependence of calcium + /// sensitivity (shifts CaRef/[Ca2+]T50 with sarcomere stretch) [-] + double beta1; }; From 36eb05c8dac486b1330336548e8141503ac6852a Mon Sep 17 00:00:00 2001 From: KB Date: Tue, 25 Aug 2026 10:43:12 -0700 Subject: [PATCH 05/13] removed unused functions and edited comments --- .../solver/active_stress_land_niederer.cpp | 14 ----- .../solver/active_stress_land_niederer.h | 56 +++++++++++++------ 2 files changed, 39 insertions(+), 31 deletions(-) diff --git a/Code/Source/solver/active_stress_land_niederer.cpp b/Code/Source/solver/active_stress_land_niederer.cpp index 440476e7e..0496f929b 100644 --- a/Code/Source/solver/active_stress_land_niederer.cpp +++ b/Code/Source/solver/active_stress_land_niederer.cpp @@ -140,18 +140,4 @@ double LandNiederer::compute_active_tension_local(const Vector &state) c return LFac * (Tref/rs) * ((ZETAS+1.0) * XS + (ZETAW) * XW); } - -double LandNiederer::compute_active_stiffness_local(const Vector &state, - const double fiber_stretch, - const double fiber_stretch_rate) const { - const double XW = std::max(0.0, state[1]); - const double XS = std::max(0.0, state[3]); - const double LFac = state[6]; - const double Aw = Aeff * rs/((1.0 - rs) * rw + rs); - const double As = Aw; - - return LFac * (Tref/rs) * (As*XS + Aw*XW); -} - - REGISTER_ACTIVE_STRESS_MODEL("LandNiederer", LandNiederer); \ No newline at end of file diff --git a/Code/Source/solver/active_stress_land_niederer.h b/Code/Source/solver/active_stress_land_niederer.h index 669bca775..433575277 100644 --- a/Code/Source/solver/active_stress_land_niederer.h +++ b/Code/Source/solver/active_stress_land_niederer.h @@ -6,6 +6,31 @@ #include "ActiveStressODE.h" +/** + * @brief Land-Niederer active stress model. + * + * This class implements the phenomenological Land-Niederer active stress model + * [1], which represents the cross-bridge cycling and calcium-troponin binding + * kinetics in a single set of ODEs. The cross-bridges are modeled as a three-state + * system, with transitions between unbound (U), pre-power stroke (W), and post-power + * stroke (S) states. The blocked state (B) represents the tropomyosin blocking + * the actin binding sites, and the unblocked state (U) repesents the availability + * of these sites for cross-bridge binding. The unblocking of the tropomyosin is + * regulated by the calcium-troponin binding state (CaTRPN) and the fiber stretch. + * The distortion of the cross-bridges are also represented for the pre-power stroke + * and post-power stroke states. The active tension is computed as a function of the + * cross-bridge states and their distortions, scaled by a reference tension Tref. + * This is given as + * @f[ + * \Tact = \frac{Tref}{r_s} \left[ (1 + \zeta_S) X_S + \zeta_W X_W \right]\;, + * @f] + * where @f$X_S@f$ and @f$X_W@f$ are the fractions of cross-bridges in the post-power + * stroke and pre-power stroke states, respectively. @f$\zeta_S@f$ and @f$\zeta_W@f$ + * are the distortions, and @f$r_s@f$ is the steady-state duty ratio of the + * cross-bridges. + * **References**: + * 1. [Land, Niederer (2017)](https://doi.org/10.1016/j.yjmcc.2017.03.008) + */ class LandNiederer : public ActiveStressODE { public: /// Model label, used for factory registration and XML selection. @@ -21,7 +46,7 @@ class LandNiederer : public ActiveStressODE { Parameters() : ActiveStressODE::Parameters(label) { constexpr bool required = true; - // Reference values: Land et al. 2017 + // Reference values calibrated in Land 2017 add_parameter("CaRef", 0.805, required); add_parameter("eta_Tm", 5.0, required); add_parameter("k_uw", 0.182, required); @@ -76,6 +101,10 @@ class LandNiederer : public ActiveStressODE { */ virtual void init_local(Vector &state) const override; + /** + * @brief Advance in time for a single node. + * + */ virtual void advance_time_step_local(const double t, const double dt, const double calcium, const double fiber_stretch, @@ -95,19 +124,11 @@ class LandNiederer : public ActiveStressODE { virtual double compute_active_tension_local(const Vector &state) const override; - /** - * @brief Compute the active stiffness for a single node. - */ - virtual double - compute_active_stiffness_local(const Vector &state, - const double fiber_stretch, - const double fiber_stretch_rate) const override; - /// @name Model parameters. /// @{ /// Reference intracellular calcium concentration giving half-maximal - /// troponin C saturation, @f$[Ca^{2+}]_{T50,ref}@f$, used (with length - /// dependence via beta1) in the CaTRPN binding ODE [uM] + /// troponin C saturation, @f$[Ca^{2+}]_{T50,ref}@f$, used in the CaTRPN + /// binding ODE [uM] double CaRef; /// Cooperativity (Hill) exponent @f$n_{Tm}@f$ for the tropomyosin @@ -115,18 +136,19 @@ class LandNiederer : public ActiveStressODE { /// activation by CaTRPN [-] double eta_Tm; - /// Rate constant for transition from unbound (U) to pre-power stroke (W) state [1/ms] + /// Rate constant for transition from unbound (U) to pre-power stroke (W) + /// state [1/ms] double k_uw; - /// Rate constant for transition from pre-power stroke (W) to post-power stroke (S) state [1/ms] + /// Rate constant for transition from pre-power stroke (W) to post-power + /// stroke (S) state [1/ms] double k_ws; /// Maximum observable tension @f$\Tref@f$ at resting length [kPa] double Tref; /// Rate constant @f$k_{TRPN}@f$ governing calcium binding/unbinding - /// kinetics of troponin C (the "pace" of CaTRPN relaxation to its - /// steady-state value) [1/ms] + /// kinetics of troponin C [1/ms] double k_TRPN; /// Cooperativity (Hill) exponent @f$n_{TRPN}@f$ of the calcium-troponin C @@ -177,8 +199,8 @@ class LandNiederer : public ActiveStressODE { /// Coefficient @f$\beta_1@f$ controlling the length-dependence of calcium /// sensitivity (shifts CaRef/[Ca2+]T50 with sarcomere stretch) [-] double beta1; -}; - + /// @} +}; #endif \ No newline at end of file From d4a627775e3a204ccd65edef147965248e35a977 Mon Sep 17 00:00:00 2001 From: KB Date: Tue, 25 Aug 2026 13:40:16 -0700 Subject: [PATCH 06/13] removed stabilization term (deferred to later issue) --- Code/Source/solver/ActiveStress.cpp | 11 +---------- Code/Source/solver/ActiveStress.h | 11 ----------- Code/Source/solver/Parameters.cpp | 9 --------- Code/Source/solver/Parameters.h | 6 ------ .../electromechanics/slab_LandNiederer/solver.xml | 1 - 5 files changed, 1 insertion(+), 37 deletions(-) diff --git a/Code/Source/solver/ActiveStress.cpp b/Code/Source/solver/ActiveStress.cpp index 32084889f..40a158802 100644 --- a/Code/Source/solver/ActiveStress.cpp +++ b/Code/Source/solver/ActiveStress.cpp @@ -14,8 +14,6 @@ void ActiveStress::read_parameters(const ActiveStressParameters ¶ms) { eta_s = params.get_eta_s(); eta_n = params.get_eta_n(); - use_stabilization = params.get_use_stabilization(); - read_model_specific_parameters( params.get_parameters(params.get_model_name())); } @@ -25,7 +23,6 @@ void ActiveStress::distribute_parameters(const CmMod &cm_mod, cm.bcast(cm_mod, &eta_f); cm.bcast(cm_mod, &eta_s); cm.bcast(cm_mod, &eta_n); - cm.bcast(cm_mod, &use_stabilization); distribute_model_specific_parameters(cm_mod, cm); } @@ -43,15 +40,9 @@ void ActiveStress::init(const unsigned int tnNo) { } active_tension.resize(tnNo); - raw_active_tension.resize(tnNo); - previous_fiber_stretch.resize(tnNo); - has_previous_fiber_stretch.resize(tnNo); for (unsigned int i = 0; i < tnNo; ++i) { - active_tension[i] = 0.0; - raw_active_tension[i] = 0.0; - previous_fiber_stretch[i] = 1.0; - has_previous_fiber_stretch[i] = 0; + active_tension[i] = 0.0; } } diff --git a/Code/Source/solver/ActiveStress.h b/Code/Source/solver/ActiveStress.h index fea5f31c4..29a419cf6 100644 --- a/Code/Source/solver/ActiveStress.h +++ b/Code/Source/solver/ActiveStress.h @@ -275,17 +275,6 @@ class ActiveStress { /// Active tension at every node. Vector active_tension; - /// Raw active tension at every node, before stabilization is applied. - Vector raw_active_tension; - - /// Previous fiber stretch at every node, used for stabilization. - Vector previous_fiber_stretch; - - /// Whether to apply stabilization to the active tension. - bool use_stabilization = false; - - Vector has_previous_fiber_stretch; - /// Active tension coefficient along the fiber direction. double eta_f; diff --git a/Code/Source/solver/Parameters.cpp b/Code/Source/solver/Parameters.cpp index 5d84c260e..f0922a258 100644 --- a/Code/Source/solver/Parameters.cpp +++ b/Code/Source/solver/Parameters.cpp @@ -1910,12 +1910,7 @@ const std::string ActiveStressParameters::xml_element_name = "Active_stress"; ActiveStressParameters::ActiveStressParameters() { model_name = Parameter("Model", "", true); - use_stabilization = - Parameter("Use_stabilization", false, false); - set_parameter("Model", "", /* required = */ true, model_name); - set_parameter("Use_stabilization", false, /* required = */ false, - use_stabilization); ActiveStressFactory::visit( [this](const std::string &name, const ActiveStress &model) { @@ -1989,10 +1984,6 @@ double ActiveStressParameters::get_eta_n() const { return directional_distribution.sheet_normal_direction.value(); } -bool ActiveStressParameters::get_use_stabilization() const { - return use_stabilization.value(); -} - const ActiveStressModelParameters & ActiveStressParameters::get_parameters(const std::string &model_name) const { if (active_stress_models.count(model_name) == 0) { diff --git a/Code/Source/solver/Parameters.h b/Code/Source/solver/Parameters.h index aa3e7a627..94d3b677b 100644 --- a/Code/Source/solver/Parameters.h +++ b/Code/Source/solver/Parameters.h @@ -1534,9 +1534,6 @@ class ActiveStressParameters : public ParameterLists { /// Get the active tension coefficient along sheet normals. double get_eta_n() const; - /// Get the stabilization flag for Regazzoni-style active tension stabilization. - bool get_use_stabilization() const; - /// Get the parameters for a given active stress model. const ActiveStressModelParameters & get_parameters(const std::string &model_name) const; @@ -1548,9 +1545,6 @@ class ActiveStressParameters : public ParameterLists { /// Parameter for the model name. Parameter model_name; - /// Flag for Regazzoni-style active tension stabilization - Parameter use_stabilization; - /// Parameters for the directional distribution of active tension. DirectionalDistributionParameters directional_distribution; diff --git a/tests/cases/electromechanics/slab_LandNiederer/solver.xml b/tests/cases/electromechanics/slab_LandNiederer/solver.xml index ccf80b824..bf7ca3985 100644 --- a/tests/cases/electromechanics/slab_LandNiederer/solver.xml +++ b/tests/cases/electromechanics/slab_LandNiederer/solver.xml @@ -134,7 +134,6 @@ LandNiederer - false 0.7 From 752feee5d5d2164096b7f750036778b9319a8634 Mon Sep 17 00:00:00 2001 From: KB Date: Tue, 25 Aug 2026 14:22:26 -0700 Subject: [PATCH 07/13] Moved Land-Niederer CI test --- tests/cases/electromechanics/slab/README.md | 21 +++++++++++++++++-- .../solver_LandNiederer.xml} | 6 +++--- .../slab_LandNiederer/mesh/X0.vtp | 3 --- .../slab_LandNiederer/mesh/X1.vtp | 3 --- .../slab_LandNiederer/mesh/volume.vtu | 3 --- 5 files changed, 22 insertions(+), 14 deletions(-) rename tests/cases/electromechanics/{slab_LandNiederer/solver.xml => slab/solver_LandNiederer.xml} (97%) delete mode 100644 tests/cases/electromechanics/slab_LandNiederer/mesh/X0.vtp delete mode 100644 tests/cases/electromechanics/slab_LandNiederer/mesh/X1.vtp delete mode 100644 tests/cases/electromechanics/slab_LandNiederer/mesh/volume.vtu diff --git a/tests/cases/electromechanics/slab/README.md b/tests/cases/electromechanics/slab/README.md index 383876bc6..99afb1f17 100755 --- a/tests/cases/electromechanics/slab/README.md +++ b/tests/cases/electromechanics/slab/README.md @@ -2,15 +2,16 @@ # **Problem Description** Simulate cardiac electromechanics on a slab of myocardial tissue. This -directory contains two solver configurations that share the same geometry and +directory contains three solver configurations that share the same geometry and electrophysiology setup but differ in the active-stress model: | Configuration file | Active-stress model | |------------------------------|---------------------| | `solver_NashPanfilov.xml` | Nash-Panfilov | | `solver_Regazzoni.xml` | RDQ20-MF (Regazzoni)| +| `solver_LandNiederer.xml` | Land-Niederer | -Both configurations couple cardiac electrophysiology (`CEP`) to solid mechanics +The configurations couple cardiac electrophysiology (`CEP`) to solid mechanics (`struct`), reproducing the geometry and stimulation setting of the Niederer electrophysiology benchmark [1] with the addition of active contraction and finite-strain mechanics. @@ -108,6 +109,18 @@ using the calcium and sarcomere-length inputs from this one-step test. The remai fields in the VTU serve as integrated svMultiPhysics regression references and were not independently validated by the RDQ20-MF reference code. +## Land-Niederer variant (`solver_LandNiederer.xml`) + +Active contraction is driven by the calcium concentration computed by the +electrophysiology model, through the Land-Niederer active-stress model [7]. The +model parameters are those calibrated to human ventricular cardiomyocyte data. + +### Validation + +The active tension fields in `result_LandNiederer_001.vtu` were validated +node-by-node against the MATLAB reference implementation provided +[here](https://www.cemrg.co.uk/models). + ## References [1] S. A. Niederer, E. Kerfoot, A. P. Benson, et al. Verification of cardiac tissue @@ -133,3 +146,7 @@ study reentrant cardiac arrhythmias. Progress in Biophysics and Molecular Biolog [6] F. Regazzoni, L. Dede', and A. Quarteroni. Biophysically detailed mathematical models of multiscale cardiac active mechanics. PLOS Computational Biology, 16(10):e1008294, 2020. + +[7] S. Land, S. Park-Holohan, N. P. Smith, et al. A model of cardiac contraction +based on novel measurements of tension development in human cardiomyocytes. +Journal of Molecular and Cellular Cardiology, 106:68-83, apr 2017. \ No newline at end of file diff --git a/tests/cases/electromechanics/slab_LandNiederer/solver.xml b/tests/cases/electromechanics/slab/solver_LandNiederer.xml similarity index 97% rename from tests/cases/electromechanics/slab_LandNiederer/solver.xml rename to tests/cases/electromechanics/slab/solver_LandNiederer.xml index bf7ca3985..26cf3862e 100644 --- a/tests/cases/electromechanics/slab_LandNiederer/solver.xml +++ b/tests/cases/electromechanics/slab/solver_LandNiederer.xml @@ -4,13 +4,13 @@ false 3 10 - 1.0 + 0.1 0.50 STOP_SIM true result - 10 + 1 0 1000 @@ -142,7 +142,7 @@ - FE + RK4 0.000805 5.0 diff --git a/tests/cases/electromechanics/slab_LandNiederer/mesh/X0.vtp b/tests/cases/electromechanics/slab_LandNiederer/mesh/X0.vtp deleted file mode 100644 index eaebcc0eb..000000000 --- a/tests/cases/electromechanics/slab_LandNiederer/mesh/X0.vtp +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:3c0889ad4a8a4659996309b5bba825c1ac5c3f03d2faca2e84e48f1e035ca240 -size 5219 diff --git a/tests/cases/electromechanics/slab_LandNiederer/mesh/X1.vtp b/tests/cases/electromechanics/slab_LandNiederer/mesh/X1.vtp deleted file mode 100644 index 07ad82435..000000000 --- a/tests/cases/electromechanics/slab_LandNiederer/mesh/X1.vtp +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:8248936f54f9537dcad8e2344beb804485d0f5cbfb960f0a8b025310e1af99ae -size 5128 diff --git a/tests/cases/electromechanics/slab_LandNiederer/mesh/volume.vtu b/tests/cases/electromechanics/slab_LandNiederer/mesh/volume.vtu deleted file mode 100644 index dd12aca89..000000000 --- a/tests/cases/electromechanics/slab_LandNiederer/mesh/volume.vtu +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:2d9c504f6d9dc221739d107e19414233cdda640d73647eba03902a1605274149 -size 421967 From 518c716c6d4e43da6d6865b7d17d8222ed79d456 Mon Sep 17 00:00:00 2001 From: KB Date: Tue, 25 Aug 2026 16:25:47 -0700 Subject: [PATCH 08/13] fixed function signature and added disable feedback parameter --- .../solver/active_stress_land_niederer.cpp | 25 +++++++++++++------ .../solver/active_stress_land_niederer.h | 18 ++++++++++--- .../slab/solver_LandNiederer.xml | 12 +++++---- 3 files changed, 40 insertions(+), 15 deletions(-) diff --git a/Code/Source/solver/active_stress_land_niederer.cpp b/Code/Source/solver/active_stress_land_niederer.cpp index 0496f929b..54f1e82f5 100644 --- a/Code/Source/solver/active_stress_land_niederer.cpp +++ b/Code/Source/solver/active_stress_land_niederer.cpp @@ -24,6 +24,12 @@ void LandNiederer::read_model_specific_parameters( Aeff = params.get_scalar("Aeff"); beta0 = params.get_scalar("beta0"); beta1 = params.get_scalar("beta1"); + + disable_force_strain_rate_feedback_ = + params.get_bool("Disable_force_strain_rate_feedback"); + + // With the feedback disabled the stretch rate computation can be skipped. + needs_fiber_stretch_rate_ = !disable_force_strain_rate_feedback_; } void LandNiederer::distribute_model_specific_parameters(const CmMod &cm_mod, @@ -47,6 +53,8 @@ void LandNiederer::distribute_model_specific_parameters(const CmMod &cm_mod, cm.bcast(cm_mod, &Aeff); cm.bcast(cm_mod, &beta0); cm.bcast(cm_mod, &beta1); + cm.bcast(cm_mod, &disable_force_strain_rate_feedback_); + cm.bcast(cm_mod, &needs_fiber_stretch_rate_); } void LandNiederer::init_local(Vector &state) const { @@ -118,19 +126,22 @@ Vector LandNiederer::getf(const double t, const Vector &state, const double gamma_su = gamma_s * std::max(term1, term2); f[3] = k_ws*XW - k_su*XS - gamma_su*XS; - // Distortion-Decay Kinetics - const double Aw = Aeff * rs/((1.0 - rs) * rw + rs); + // Distortion-Decay Kinetics + const double stretch_rate = + disable_force_strain_rate_feedback_ ? 0.0 : fiber_stretch_rate; + const double Aw = Aeff * rs/((1.0 - rs) * rw + rs); const double cw = phi * k_uw * (1.0 - rs)*(1.0 - rw) / ((1.0 - rs)*rw); - f[4] = Aw * fiber_stretch_rate - cw*ZETAW; + f[4] = Aw * stretch_rate - cw*ZETAW; - const double As = Aw; - const double cs = phi * k_ws * (1.0 - rs)*rw/rs; - f[5] = As * fiber_stretch_rate - cs*ZETAS; + const double As = Aw; + const double cs = phi * k_ws * (1.0 - rs)*rw/rs; + f[5] = As * stretch_rate - cs*ZETAS; return f; } -double LandNiederer::compute_active_tension_local(const Vector &state) const { +double LandNiederer::compute_active_tension_local(const Vector &state, + const double fiber_stretch) const { const double XW = std::max(0.0, state[1]); const double XS = std::max(0.0, state[3]); const double ZETAW = state[4]; diff --git a/Code/Source/solver/active_stress_land_niederer.h b/Code/Source/solver/active_stress_land_niederer.h index 433575277..8a462e51c 100644 --- a/Code/Source/solver/active_stress_land_niederer.h +++ b/Code/Source/solver/active_stress_land_niederer.h @@ -63,14 +63,18 @@ class LandNiederer : public ActiveStressODE { add_parameter("phi", 2.23, required); add_parameter("Aeff", 25.0, required); add_parameter("beta0", 2.3, required); - add_parameter("beta1", -2.4, required); + add_parameter("beta1", -2.4, required); + + add_parameter("Disable_force_strain_rate_feedback", false, !required); } }; /** * @brief Constructor. */ - LandNiederer() : ActiveStressODE(7) {} + LandNiederer() : ActiveStressODE(/* n_state_variables = */ 7, + /* needs_fiber_stretch = */ true, + /* needs_fiber_stretch_rate = */ true) {} /** * @brief Construct an instance of model parameters. @@ -122,7 +126,8 @@ class LandNiederer : public ActiveStressODE { * @brief Compute the active tension for a single node. */ virtual double - compute_active_tension_local(const Vector &state) const override; + compute_active_tension_local(const Vector &state, + const double fiber_stretch) const override; /// @name Model parameters. /// @{ @@ -200,6 +205,13 @@ class LandNiederer : public ActiveStressODE { /// sensitivity (shifts CaRef/[Ca2+]T50 with sarcomere stretch) [-] double beta1; + /// Controls force-strain-rate feedback in the cross-bridge distortion ODEs. + /// - @c false (default): the ZETAW/ZETAS distortion states respond to the + /// fiber stretch rate as in Land, Niederer (2017). + /// - @c true: the fiber stretch rate contribution is zeroed out, disabling + /// the feedback. + bool disable_force_strain_rate_feedback_ = false; + /// @} }; diff --git a/tests/cases/electromechanics/slab/solver_LandNiederer.xml b/tests/cases/electromechanics/slab/solver_LandNiederer.xml index 26cf3862e..a7ea01f58 100644 --- a/tests/cases/electromechanics/slab/solver_LandNiederer.xml +++ b/tests/cases/electromechanics/slab/solver_LandNiederer.xml @@ -4,7 +4,7 @@ false 3 10 - 0.1 + 1.0 0.50 STOP_SIM @@ -96,7 +96,7 @@ true - + fsils @@ -159,8 +159,10 @@ 0.615 2.23 25.0 - 0.0 - 0.0 + 2.3 + -2.4 + + true @@ -184,7 +186,7 @@ fsils 1e-12 - 1e-14 + 1e-12 1000 From 03b671b5f1e51cf2094e5dbe7f3f4ae705062547 Mon Sep 17 00:00:00 2001 From: KB Date: Thu, 3 Sep 2026 11:35:01 -0400 Subject: [PATCH 09/13] Added CI Test --- tests/cases/electromechanics/slab/README.md | 33 ++++++++++++++----- .../slab/result_LandNiederer_001.vtu | 3 ++ .../slab/solver_LandNiederer.xml | 6 ++-- tests/test_electromechanics.py | 2 +- 4 files changed, 32 insertions(+), 12 deletions(-) create mode 100644 tests/cases/electromechanics/slab/result_LandNiederer_001.vtu diff --git a/tests/cases/electromechanics/slab/README.md b/tests/cases/electromechanics/slab/README.md index 99afb1f17..d5167435a 100755 --- a/tests/cases/electromechanics/slab/README.md +++ b/tests/cases/electromechanics/slab/README.md @@ -88,6 +88,28 @@ are not prescribed by the RDQ20-MF model itself. **Regression reference:** `result_Regazzoni_001.vtu` +## Land-Niederer variant (`solver_LandNiederer.xml`) + +Active contraction is driven by the calcium concentration computed by the +electrophysiology model, through the Land-Niederer active-stress model [7]. The +model parameters are those calibrated to human ventricular cardiomyocyte data. The +scalar active tension is distributed along the fiber, sheet, and sheet-normal +directions using the same directional weights as the Nash-Panfilov variant. + +``` + + LandNiederer + + 0.7 + 0.2 + 0.1 + + ... + +``` + +**Regression reference:** `result_LandNiederer_001.vtu` + ### Validation svMultiPhysics stores `T_act = a_XB * (μ_P^1 + μ_N^1) * φ(SL)` — the scalar @@ -109,17 +131,10 @@ using the calcium and sarcomere-length inputs from this one-step test. The remai fields in the VTU serve as integrated svMultiPhysics regression references and were not independently validated by the RDQ20-MF reference code. -## Land-Niederer variant (`solver_LandNiederer.xml`) - -Active contraction is driven by the calcium concentration computed by the -electrophysiology model, through the Land-Niederer active-stress model [7]. The -model parameters are those calibrated to human ventricular cardiomyocyte data. - -### Validation - The active tension fields in `result_LandNiederer_001.vtu` were validated node-by-node against the MATLAB reference implementation provided -[here](https://www.cemrg.co.uk/models). +[here](https://www.cemrg.co.uk/models). A custom script with a python wrapper +for the MATLAB function is available in this [repository](https://github.com/kko27/generate_ref_solution.git). ## References diff --git a/tests/cases/electromechanics/slab/result_LandNiederer_001.vtu b/tests/cases/electromechanics/slab/result_LandNiederer_001.vtu new file mode 100644 index 000000000..1af9d63e5 --- /dev/null +++ b/tests/cases/electromechanics/slab/result_LandNiederer_001.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:804a6ecf2b5bacb1a8e55392a6302138652e58f883a17ee10c99e4397cad334d +size 1470461 diff --git a/tests/cases/electromechanics/slab/solver_LandNiederer.xml b/tests/cases/electromechanics/slab/solver_LandNiederer.xml index a7ea01f58..89f3f80fc 100644 --- a/tests/cases/electromechanics/slab/solver_LandNiederer.xml +++ b/tests/cases/electromechanics/slab/solver_LandNiederer.xml @@ -3,7 +3,7 @@ false 3 - 10 + 1 1.0 0.50 STOP_SIM @@ -179,6 +179,8 @@ true true true + + true @@ -186,7 +188,7 @@ fsils 1e-12 - 1e-12 + 1e-14 1000 diff --git a/tests/test_electromechanics.py b/tests/test_electromechanics.py index e4f63bbd0..6f79c1749 100644 --- a/tests/test_electromechanics.py +++ b/tests/test_electromechanics.py @@ -23,7 +23,7 @@ ] -@pytest.mark.parametrize("model", ["NashPanfilov", "Regazzoni"]) +@pytest.mark.parametrize("model", ["NashPanfilov", "Regazzoni", "LandNiederer"]) def test_slab(model, n_proc): run_with_reference(base_folder, "slab", fields, n_proc, t_max=1, name_inp=f"solver_{model}.xml", From 0d22ab552945a323448fd0b2d66a4d42f52429db Mon Sep 17 00:00:00 2001 From: KB Date: Thu, 3 Sep 2026 12:06:33 -0400 Subject: [PATCH 10/13] Rename Land-Niederer active stress model to CamelCase convention Follows the file/class naming introduced upstream (ActiveStressNashPanfilov, ActiveStressRegazzoni, ...): - active_stress_land_niederer.{cpp,h} -> ActiveStressLandNiederer.{cpp,h} - class LandNiederer -> ActiveStressLandNiederer The XML model label ("LandNiederer") is unchanged. Co-Authored-By: Claude Sonnet 5 --- ...niederer.cpp => ActiveStressLandNiederer.cpp} | 16 ++++++++-------- ...and_niederer.h => ActiveStressLandNiederer.h} | 4 ++-- Code/Source/solver/CMakeLists.txt | 2 +- .../test_electromechanics_models.cpp | 4 ++-- 4 files changed, 13 insertions(+), 13 deletions(-) rename Code/Source/solver/{active_stress_land_niederer.cpp => ActiveStressLandNiederer.cpp} (89%) rename Code/Source/solver/{active_stress_land_niederer.h => ActiveStressLandNiederer.h} (98%) diff --git a/Code/Source/solver/active_stress_land_niederer.cpp b/Code/Source/solver/ActiveStressLandNiederer.cpp similarity index 89% rename from Code/Source/solver/active_stress_land_niederer.cpp rename to Code/Source/solver/ActiveStressLandNiederer.cpp index 54f1e82f5..0834db03e 100644 --- a/Code/Source/solver/active_stress_land_niederer.cpp +++ b/Code/Source/solver/ActiveStressLandNiederer.cpp @@ -1,9 +1,9 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the // University of California, and others. SPDX-License-Identifier: BSD-3-Clause -#include "active_stress_land_niederer.h" +#include "ActiveStressLandNiederer.h" -void LandNiederer::read_model_specific_parameters( +void ActiveStressLandNiederer::read_model_specific_parameters( const ActiveStressModelParameters ¶ms) { ActiveStressODE::read_model_specific_parameters(params); @@ -32,7 +32,7 @@ void LandNiederer::read_model_specific_parameters( needs_fiber_stretch_rate_ = !disable_force_strain_rate_feedback_; } -void LandNiederer::distribute_model_specific_parameters(const CmMod &cm_mod, +void ActiveStressLandNiederer::distribute_model_specific_parameters(const CmMod &cm_mod, const cmType &cm) { ActiveStressODE::distribute_model_specific_parameters(cm_mod, cm); @@ -57,7 +57,7 @@ void LandNiederer::distribute_model_specific_parameters(const CmMod &cm_mod, cm.bcast(cm_mod, &needs_fiber_stretch_rate_); } -void LandNiederer::init_local(Vector &state) const { +void ActiveStressLandNiederer::init_local(Vector &state) const { state[0] = 1.0; state[1] = 0.0; state[2] = 0.0; @@ -68,7 +68,7 @@ void LandNiederer::init_local(Vector &state) const { state[6] = std::max(0.0, 1.0 + beta0 * (lambda + std::min(0.87, lambda) - 1.87)); } -void LandNiederer::advance_time_step_local(const double t, const double dt, +void ActiveStressLandNiederer::advance_time_step_local(const double t, const double dt, const double calcium, const double fiber_stretch, const double fiber_stretch_rate, @@ -90,7 +90,7 @@ void LandNiederer::advance_time_step_local(const double t, const double dt, state[6] = std::max(0.0, 1.0 + beta0 * (lambda + std::min(0.87, lambda) - 1.87)); } -Vector LandNiederer::getf(const double t, const Vector &state, +Vector ActiveStressLandNiederer::getf(const double t, const Vector &state, const double calcium, const double fiber_stretch, const double fiber_stretch_rate) const { @@ -140,7 +140,7 @@ Vector LandNiederer::getf(const double t, const Vector &state, return f; } -double LandNiederer::compute_active_tension_local(const Vector &state, +double ActiveStressLandNiederer::compute_active_tension_local(const Vector &state, const double fiber_stretch) const { const double XW = std::max(0.0, state[1]); const double XS = std::max(0.0, state[3]); @@ -151,4 +151,4 @@ double LandNiederer::compute_active_tension_local(const Vector &state, return LFac * (Tref/rs) * ((ZETAS+1.0) * XS + (ZETAW) * XW); } -REGISTER_ACTIVE_STRESS_MODEL("LandNiederer", LandNiederer); \ No newline at end of file +REGISTER_ACTIVE_STRESS_MODEL("LandNiederer", ActiveStressLandNiederer); \ No newline at end of file diff --git a/Code/Source/solver/active_stress_land_niederer.h b/Code/Source/solver/ActiveStressLandNiederer.h similarity index 98% rename from Code/Source/solver/active_stress_land_niederer.h rename to Code/Source/solver/ActiveStressLandNiederer.h index 8a462e51c..46bf80530 100644 --- a/Code/Source/solver/active_stress_land_niederer.h +++ b/Code/Source/solver/ActiveStressLandNiederer.h @@ -31,7 +31,7 @@ * **References**: * 1. [Land, Niederer (2017)](https://doi.org/10.1016/j.yjmcc.2017.03.008) */ -class LandNiederer : public ActiveStressODE { +class ActiveStressLandNiederer : public ActiveStressODE { public: /// Model label, used for factory registration and XML selection. static inline const std::string label = "LandNiederer"; @@ -72,7 +72,7 @@ class LandNiederer : public ActiveStressODE { /** * @brief Constructor. */ - LandNiederer() : ActiveStressODE(/* n_state_variables = */ 7, + ActiveStressLandNiederer() : ActiveStressODE(/* n_state_variables = */ 7, /* needs_fiber_stretch = */ true, /* needs_fiber_stretch_rate = */ true) {} diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt index 7f8d1d775..fe423e708 100644 --- a/Code/Source/solver/CMakeLists.txt +++ b/Code/Source/solver/CMakeLists.txt @@ -269,7 +269,7 @@ set(CSRCS ActiveStressODE.cpp ActiveStressNashPanfilov.cpp ActiveStressRegazzoni.cpp - active_stress_land_niederer.cpp + ActiveStressLandNiederer.cpp SPLIT.c diff --git a/tests/unitTests/ionic_model_tests/test_electromechanics_models.cpp b/tests/unitTests/ionic_model_tests/test_electromechanics_models.cpp index 941e34a70..eea548704 100644 --- a/tests/unitTests/ionic_model_tests/test_electromechanics_models.cpp +++ b/tests/unitTests/ionic_model_tests/test_electromechanics_models.cpp @@ -30,9 +30,9 @@ */ #include "../test_common.h" -#include "active_stress_land_niederer.h" +#include "ActiveStressLandNiederer.h" -class TestableLandNiederer : public LandNiederer { +class TestableLandNiederer : public ActiveStressLandNiederer { public: Vector evaluate_getf(const double t, const Vector &state, const double calcium, From 40fe39599e007abbf32862dfa4c29fae63d98d5a Mon Sep 17 00:00:00 2001 From: KB Date: Tue, 22 Sep 2026 13:16:05 -0700 Subject: [PATCH 11/13] added Ca scaling factor and updated .xml file --- Code/Source/solver/ActiveStressLandNiederer.cpp | 4 +++- Code/Source/solver/ActiveStressLandNiederer.h | 5 +++++ tests/cases/electromechanics/slab/solver_LandNiederer.xml | 8 +++++--- 3 files changed, 13 insertions(+), 4 deletions(-) diff --git a/Code/Source/solver/ActiveStressLandNiederer.cpp b/Code/Source/solver/ActiveStressLandNiederer.cpp index 0834db03e..b4295df4b 100644 --- a/Code/Source/solver/ActiveStressLandNiederer.cpp +++ b/Code/Source/solver/ActiveStressLandNiederer.cpp @@ -7,6 +7,7 @@ void ActiveStressLandNiederer::read_model_specific_parameters( const ActiveStressModelParameters ¶ms) { ActiveStressODE::read_model_specific_parameters(params); + calcium_scaling_factor = params.get_scalar("calcium_scaling_factor"); CaRef = params.get_scalar("CaRef"); eta_Tm = params.get_scalar("eta_Tm"); k_uw = params.get_scalar("k_uw"); @@ -36,6 +37,7 @@ void ActiveStressLandNiederer::distribute_model_specific_parameters(const CmMod const cmType &cm) { ActiveStressODE::distribute_model_specific_parameters(cm_mod, cm); + cm.bcast(cm_mod, &calcium_scaling_factor); cm.bcast(cm_mod, &CaRef); cm.bcast(cm_mod, &eta_Tm); cm.bcast(cm_mod, &k_uw); @@ -117,7 +119,7 @@ Vector ActiveStressLandNiederer::getf(const double t, const Vector RK4 - 0.000805 + + 1000.0 + 0.805 5.0 0.182 0.012 @@ -157,10 +159,10 @@ 0.25 0.0085 0.615 - 2.23 + 2.23 25.0 2.3 - -2.4 + -2.4 true From e735ce3d9d505af32c66e72c9e60e482b04a4f83 Mon Sep 17 00:00:00 2001 From: KB Date: Wed, 23 Sep 2026 13:52:14 -0700 Subject: [PATCH 12/13] added unit test case for LN Model --- .../Testing/Temporary/CTestCostData.txt | 1 - .../unitTests/Testing/Temporary/LastTest.log | 3 - .../test_active_stress_land_niederer.cpp | 38 +++++++ .../test_electromechanics_models.cpp | 104 ------------------ .../active_stress_land_niederer_twitch.csv | 12 ++ .../active_stress/land_niederer/README.md | 24 ++++ 6 files changed, 74 insertions(+), 108 deletions(-) delete mode 100644 tests/unitTests/Testing/Temporary/CTestCostData.txt delete mode 100644 tests/unitTests/Testing/Temporary/LastTest.log create mode 100644 tests/unitTests/active_stress_tests/test_active_stress_land_niederer.cpp delete mode 100644 tests/unitTests/ionic_model_tests/test_electromechanics_models.cpp create mode 100644 tests/unitTests/reference_data/active_stress_land_niederer_twitch.csv create mode 100644 tests/unitTests/reference_generators/active_stress/land_niederer/README.md diff --git a/tests/unitTests/Testing/Temporary/CTestCostData.txt b/tests/unitTests/Testing/Temporary/CTestCostData.txt deleted file mode 100644 index ed97d539c..000000000 --- a/tests/unitTests/Testing/Temporary/CTestCostData.txt +++ /dev/null @@ -1 +0,0 @@ ---- diff --git a/tests/unitTests/Testing/Temporary/LastTest.log b/tests/unitTests/Testing/Temporary/LastTest.log deleted file mode 100644 index d25f1e60c..000000000 --- a/tests/unitTests/Testing/Temporary/LastTest.log +++ /dev/null @@ -1,3 +0,0 @@ -Start testing: Jul 27 11:37 PDT ----------------------------------------------------------- -End testing: Jul 27 11:37 PDT diff --git a/tests/unitTests/active_stress_tests/test_active_stress_land_niederer.cpp b/tests/unitTests/active_stress_tests/test_active_stress_land_niederer.cpp new file mode 100644 index 000000000..a4cc64e3e --- /dev/null +++ b/tests/unitTests/active_stress_tests/test_active_stress_land_niederer.cpp @@ -0,0 +1,38 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause + +/// @file + +#include "ActiveStressLandNiederer.h" +#include "active_stress_test_helpers.h" +#include "gtest/gtest.h" + +/** + * @test Run a standalone ActiveStressLandNiederer twitch for 600 updates with + * @f$\Delta t=1\,\mathrm{ms}@f$, the prescribed calcium and fiber-stretch + * trajectories, and the multiscale Land 2017 calibration configured below + * (svMP's model defaults). + * + * The trusted reference is generated by driving the original MATLAB + * @c Land.m implementation (adapted from https://www.cemrg.co.uk/models) + * through the MATLAB engine with a manual Forward-Euler integrator that + * mirrors @c ActiveStressODE's update rule; it is not generated by svMultiPhysics. + * See this [repo](https://github.com/kko27/generate_ref_solution.git) to access the + * reference solution generator and model. See @ref ActiveStressLandNiederer + * for the model equations and state ordering: states 0-5 are the Land 2017 ODE states + * (XB, XW, CaTRPN, XS, ZETAW, ZETAS) and state 6 is the deterministic Frank-Starling + * length-dependence factor LFac. + */ +TEST(ActiveStressTrajectory, LandNiederer) { + ActiveStressLandNiederer::Parameters params; + + ActiveStressTrajectoryConfiguration configuration; + configuration.final_time = 600.0; + configuration.time_step = 1.0; + configuration.reference_csv_filename = + "active_stress_land_niederer_twitch.csv"; + + ActiveStressTrajectoryTest trajectory( + params, configuration); + trajectory.run(); +} diff --git a/tests/unitTests/ionic_model_tests/test_electromechanics_models.cpp b/tests/unitTests/ionic_model_tests/test_electromechanics_models.cpp deleted file mode 100644 index eea548704..000000000 --- a/tests/unitTests/ionic_model_tests/test_electromechanics_models.cpp +++ /dev/null @@ -1,104 +0,0 @@ -/* Copyright (c) Stanford University, The Regents of the University of - * California, and others. - * - * All Rights Reserved. - * - * See Copyright-SimVascular.txt for additional details. - * - * Permission is hereby granted, free of charge, to any person obtaining - * a copy of this software and associated documentation files (the - * "Software"), to deal in the Software without restriction, including - * without limitation the rights to use, copy, modify, merge, publish, - * distribute, sublicense, and/or sell copies of the Software, and to - * permit persons to whom the Software is furnished to do so, subject - * to the following conditions: - * - * The above copyright notice and this permission notice shall be included - * in all copies or substantial portions of the Software. - * - * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS - * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED - * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A - * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER - * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, - * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, - * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR - * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF - * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING - * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS - * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - */ - -#include "../test_common.h" -#include "ActiveStressLandNiederer.h" - -class TestableLandNiederer : public ActiveStressLandNiederer { -public: - Vector evaluate_getf(const double t, const Vector &state, - const double calcium, - const double fiber_stretch, - const double fiber_stretch_rate) const { - return getf(t, state, calcium, fiber_stretch, fiber_stretch_rate); - } - - void set_parameters_for_test() { - CaRef = 0.805; - eta_Tm = 5.0; - k_uw = 0.182; - k_ws = 0.012; - Tref = 120.0; - k_TRPN = 0.1; - eta_TRPN = 2.0; - k_u = 1.0; - TRPN50 = 0.35; - rw = 0.5; - rs = 0.25; - gamma_s = 0.0085; - gamma_w = 0.615; - phi = 2.23; - Aeff = 25.0; - beta0 = 2.3; - beta1 = -2.4; - } -}; - -class ElectromechanicsModelTest : public ::testing::Test { -protected: - void SetUp() override { model.set_parameters_for_test(); } - - void TearDown() override {} - - static Vector initial_state() { - return {0.8900, 0.0444, 0.2609, 0.0052, 0.0000, 0.0000}; - } - - static Vector expected_rhs() { - return {-0.0066, 0.0029, 0.0029, 0.0004, 1.0, 1.0}; - } - - TestableLandNiederer model; -}; - -TEST_F(ElectromechanicsModelTest, LandNiedererGetf) { - // Evaluate the Land-Niederer ODE right-hand side at one state and compare - // against reference values computed from the MATLAB implementation. - - const double t = 0.0; - const double calcium = 0.5040; - const double fiber_stretch = 1.0; - const double fiber_stretch_rate = 0.1; - - const Vector state = initial_state(); - const Vector expected = expected_rhs(); - - const Vector rhs = - model.evaluate_getf(t, state, calcium, fiber_stretch, - fiber_stretch_rate); - - ASSERT_EQ(rhs.size(), expected.size()); - - for (int i = 0; i < rhs.size(); ++i) { - ASSERT_NEAR(rhs[i], expected[i], 1e-3) - << "RHS mismatch at component " << i; - } -} \ No newline at end of file diff --git a/tests/unitTests/reference_data/active_stress_land_niederer_twitch.csv b/tests/unitTests/reference_data/active_stress_land_niederer_twitch.csv new file mode 100644 index 000000000..ac7a17f61 --- /dev/null +++ b/tests/unitTests/reference_data/active_stress_land_niederer_twitch.csv @@ -0,0 +1,12 @@ +step,Ta,s0,s1,s2,s3,s4,s5,s6 +0,0.0000000000000000e+00,1.0000000000000000e+00,0.0000000000000000e+00,1.5431503414220131e-03,0.0000000000000000e+00,0.0000000000000000e+00,0.0000000000000000e+00,1.0000000000000000e+00 +10,3.8835726655655308e-05,9.9999002345980692e-01,2.9941229181942714e-06,1.0517193396842484e-02,8.0907763865948556e-08,0.0000000000000000e+00,0.0000000000000000e+00,1.0000000000000000e+00 +30,2.0300434401505334e+00,5.9239113723449799e-01,1.2882747028414332e-01,4.5983992372518667e-01,4.2360968011456853e-03,-5.1401253665650515e-05,-5.1401253665650515e-05,1.0000000000000000e+00 +60,4.5116286702089390e+01,1.3635660242474829e-01,3.8474842126758285e-01,5.0941896089957439e-01,1.0076775625847574e-01,-6.6588419265878676e-03,-3.2292507987109829e-02,9.8989518395093545e-01 +99,6.0941402944749235e+01,5.6329632662556939e-01,1.6508974083429478e-01,2.8363900072384202e-01,1.4600652272940820e-01,-9.4503553504148817e-03,-8.1107979806100003e-02,9.5744613494697117e-01 +149,3.1250050187602852e+01,9.2019661964820609e-01,3.7783373861371913e-03,9.2371464543242410e-02,7.3204243866698743e-02,-4.9622131664648879e-04,-4.4722231624103201e-02,9.3101182228834278e-01 +199,1.3831174797835914e+01,9.6871352412464840e-01,6.5656925040025985e-04,3.5975407531667690e-02,3.0054736574624694e-02,3.9747232336301581e-03,1.9073206439169550e-02,9.4072464272716960e-01 +299,2.4009864451104130e+00,9.9499349091655198e-01,9.9584290598387748e-05,1.7437058100832695e-02,4.8139270752370631e-03,4.2277664579140582e-03,5.0009060755114958e-02,9.8950899148243476e-01 +349,9.4836453242079399e-01,9.9796063893219689e-01,5.2352355165063066e-05,1.6198462451437204e-02,1.9377442578147675e-03,1.7892857982496535e-04,1.9617772812423967e-02,9.9999574382061696e-01 +449,1.5914313124945889e-01,9.9961814698125151e-01,2.5461006769429786e-05,1.5334951073125025e-02,3.3144008551800805e-04,4.3808552516745440e-27,3.2616629597560521e-04,1.0000000000000000e+00 +599,1.6805285673456662e-02,9.9992369564595407e-01,2.0663311129012729e-05,1.5203845072689088e-02,3.5010987338793960e-05,5.3073494937017952e-61,6.9923499097932138e-07,1.0000000000000000e+00 diff --git a/tests/unitTests/reference_generators/active_stress/land_niederer/README.md b/tests/unitTests/reference_generators/active_stress/land_niederer/README.md new file mode 100644 index 000000000..32939fa77 --- /dev/null +++ b/tests/unitTests/reference_generators/active_stress/land_niederer/README.md @@ -0,0 +1,24 @@ +# Land-Niederer active-stress reference + +The committed reference trajectory was produced with KB Ko's [`generate_ref_solution`](https://github.com/kko27/generate_ref_solution.git) using the Land-Niederer +human ventricular cardiomyocyte parameter set for multiscale simulations. It +was not generated by svMultiPhysics. + +The reference uses 600 outer updates with `dt=1 ms`. At the start of each +update, calcium follows the normalized double-exponential transient prescribed +by `ActiveStressTrajectoryTest`: baseline `1e-4 mM`, peak `9e-4 mM`, rise and +decay constants `20 ms` and `50 ms`, and onset `10 ms`. Sarcomere length follows +the same raised-cosine trajectory, shortening from `2.2 um` at `30 ms` to +`2.134 um` at `150 ms` and recovering to `2.2 um` at `350 ms`. Its rate is the +forward difference `[SL(t+dt)-SL(t)]/dt`. + +For each update, the authors' `solve_time_step` receives time in seconds, +calcium in micromolar, sarcomere length in micrometers, and its rate in +micrometers per second. The authors' code advances the states with Forward-Euler +time steps. The Active tension is divided by 1000 to convert kPa to MPa. + +The committed CSV stores `step,Ta,s0,...,s6` at update indices +`0,10,30,60,99,149,199,299,349,449,599`, where index `N` is the result of the +update starting at `t=N ms`. No ready-to-run generator is included because +reproducing the data depends on installation of MATLAB. The trusted data remain in +`tests/unitTests/reference_data/active_stress_land_niederer_twitch.csv`. From 5ecc006e6790f29baa6d36b78fa80b67bbea32e0 Mon Sep 17 00:00:00 2001 From: KB Date: Wed, 23 Sep 2026 20:57:23 -0700 Subject: [PATCH 13/13] removed initialization redundancy --- Code/Source/solver/ActiveStress.cpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/Code/Source/solver/ActiveStress.cpp b/Code/Source/solver/ActiveStress.cpp index 40a158802..7db495dae 100644 --- a/Code/Source/solver/ActiveStress.cpp +++ b/Code/Source/solver/ActiveStress.cpp @@ -40,10 +40,6 @@ void ActiveStress::init(const unsigned int tnNo) { } active_tension.resize(tnNo); - - for (unsigned int i = 0; i < tnNo; ++i) { - active_tension[i] = 0.0; - } } void ActiveStress::advance_time_step(const double t, const double dt,