From 06557e86934eec4ac2669bfff57a86dc8f337523 Mon Sep 17 00:00:00 2001 From: Sami Bismar Date: Thu, 23 Jul 2026 12:15:39 -0500 Subject: [PATCH 01/17] Extend active-stress tension interface for fiber stretch --- Code/Source/solver/active_stress.cpp | 2 +- Code/Source/solver/active_stress.h | 6 +++++- Code/Source/solver/active_stress_nash_panfilov.cpp | 3 ++- Code/Source/solver/active_stress_nash_panfilov.h | 3 ++- Code/Source/solver/active_stress_uniform_steady.h | 3 ++- Code/Source/solver/active_stress_uniform_unsteady.cpp | 2 +- Code/Source/solver/active_stress_uniform_unsteady.h | 3 ++- 7 files changed, 15 insertions(+), 7 deletions(-) diff --git a/Code/Source/solver/active_stress.cpp b/Code/Source/solver/active_stress.cpp index d8630565d..c3ee67172 100644 --- a/Code/Source/solver/active_stress.cpp +++ b/Code/Source/solver/active_stress.cpp @@ -54,6 +54,6 @@ void ActiveStress::advance_time_step(const double t, const double dt, fiber_stretch_rate[i], state_loc); states.set_col(i, state_loc); - active_tension[i] = compute_active_tension_local(state_loc); + active_tension[i] = compute_active_tension_local(state_loc, fiber_stretch[i]); } } \ No newline at end of file diff --git a/Code/Source/solver/active_stress.h b/Code/Source/solver/active_stress.h index da5f38c26..a4427c606 100644 --- a/Code/Source/solver/active_stress.h +++ b/Code/Source/solver/active_stress.h @@ -225,9 +225,13 @@ class ActiveStress { /** * @brief Compute the active tension for a single node. + * + * @param[in] state State vector for a single node. + * @param[in] fiber_stretch Fiber stretch at the current node. */ virtual double - compute_active_tension_local(const Vector &state) const = 0; + compute_active_tension_local(const Vector &state, + const double fiber_stretch) const = 0; /// Current time. Updated whenever calling @ref advance_time_step. double time; diff --git a/Code/Source/solver/active_stress_nash_panfilov.cpp b/Code/Source/solver/active_stress_nash_panfilov.cpp index 3e23f21b2..31095da3e 100644 --- a/Code/Source/solver/active_stress_nash_panfilov.cpp +++ b/Code/Source/solver/active_stress_nash_panfilov.cpp @@ -45,7 +45,8 @@ Vector NashPanfilov::getf(const double t, const Vector &state, } double -NashPanfilov::compute_active_tension_local(const Vector &state) const { +NashPanfilov::compute_active_tension_local(const Vector &state, + const double fiber_stretch) const { return state[0]; } diff --git a/Code/Source/solver/active_stress_nash_panfilov.h b/Code/Source/solver/active_stress_nash_panfilov.h index a164be7c6..b01d6f120 100644 --- a/Code/Source/solver/active_stress_nash_panfilov.h +++ b/Code/Source/solver/active_stress_nash_panfilov.h @@ -101,7 +101,8 @@ class NashPanfilov : 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. /// @{ diff --git a/Code/Source/solver/active_stress_uniform_steady.h b/Code/Source/solver/active_stress_uniform_steady.h index 93676f9b0..ce5e2fd41 100644 --- a/Code/Source/solver/active_stress_uniform_steady.h +++ b/Code/Source/solver/active_stress_uniform_steady.h @@ -78,7 +78,8 @@ class UniformSteadyActiveStress : public ActiveStress { * @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 { return value; } diff --git a/Code/Source/solver/active_stress_uniform_unsteady.cpp b/Code/Source/solver/active_stress_uniform_unsteady.cpp index 5d02d76b2..708424085 100644 --- a/Code/Source/solver/active_stress_uniform_unsteady.cpp +++ b/Code/Source/solver/active_stress_uniform_unsteady.cpp @@ -26,7 +26,7 @@ void UniformUnsteadyActiveStress::distribute_model_specific_parameters( } double UniformUnsteadyActiveStress::compute_active_tension_local( - const Vector &state) const { + const Vector &state, const double fiber_stretch) const { return fourier_interpolation.value(time)[0]; } diff --git a/Code/Source/solver/active_stress_uniform_unsteady.h b/Code/Source/solver/active_stress_uniform_unsteady.h index 696be51ab..f28b05165 100644 --- a/Code/Source/solver/active_stress_uniform_unsteady.h +++ b/Code/Source/solver/active_stress_uniform_unsteady.h @@ -90,7 +90,8 @@ class UniformUnsteadyActiveStress : public ActiveStress { * @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; /// Toggle between ramp or Fourier transform. bool ramp; From 88702ff9a7b1659a7b016339310049f0070c7881 Mon Sep 17 00:00:00 2001 From: Sami Bismar Date: Thu, 23 Jul 2026 15:47:39 -0500 Subject: [PATCH 02/17] Add RDQ20-MF active-stress model foundation --- Code/Source/solver/CMakeLists.txt | 3 +- Code/Source/solver/active_stress_rdq20_mf.cpp | 45 +++++++ Code/Source/solver/active_stress_rdq20_mf.h | 123 ++++++++++++++++++ 3 files changed, 170 insertions(+), 1 deletion(-) create mode 100644 Code/Source/solver/active_stress_rdq20_mf.cpp create mode 100644 Code/Source/solver/active_stress_rdq20_mf.h diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt index c5ab81146..f4688f25d 100644 --- a/Code/Source/solver/CMakeLists.txt +++ b/Code/Source/solver/CMakeLists.txt @@ -234,7 +234,8 @@ set(CSRCS active_stress_uniform_unsteady.cpp active_stress_ode.cpp active_stress_nash_panfilov.cpp - + active_stress_rdq20_mf.cpp + SPLIT.c svZeroD_interface/LPNSolverInterface.h svZeroD_interface/LPNSolverInterface.cpp diff --git a/Code/Source/solver/active_stress_rdq20_mf.cpp b/Code/Source/solver/active_stress_rdq20_mf.cpp new file mode 100644 index 000000000..c5cbb9155 --- /dev/null +++ b/Code/Source/solver/active_stress_rdq20_mf.cpp @@ -0,0 +1,45 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause + +#include "active_stress_rdq20_mf.h" + +#include "Core/Exception.h" + +void RDQ20MF::read_model_specific_parameters( + const ActiveStressModelParameters ¶ms) { + // TODO: read model-specific parameters (added with the RU/XB dynamics in + // later increments). No parameters to read yet. +} + +void RDQ20MF::distribute_model_specific_parameters(const CmMod &cm_mod, + const cmType &cm) { + // TODO: distribute model-specific parameters (added with the RU/XB dynamics + // in later increments). No parameters to distribute yet. +} + +void RDQ20MF::init_local(Vector &state) const { + for (unsigned int i = 0; i < n_states; ++i) + state[i] = 0.0; + + state[ru_index(0, 0, 0, 0)] = 1.0; // == state[0] +} + +void RDQ20MF::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 { + svmp::not_implemented( + "RDQ20-MF active stress dynamics are not implemented yet: the " + "regulatory-unit and crossbridge state updates are added in a later " + "increment."); +} + +double RDQ20MF::compute_active_tension_local(const Vector &state, + const double fiber_stretch) const { + // TEMPORARY: returns zero until the active-tension formula is implemented in + // a later increment. + return 0.0; +} + +REGISTER_ACTIVE_STRESS_MODEL("RDQ20-MF", RDQ20MF); diff --git a/Code/Source/solver/active_stress_rdq20_mf.h b/Code/Source/solver/active_stress_rdq20_mf.h new file mode 100644 index 000000000..e4f557baf --- /dev/null +++ b/Code/Source/solver/active_stress_rdq20_mf.h @@ -0,0 +1,123 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause + +#ifndef ACTIVE_STRESS_RDQ20_MF_H +#define ACTIVE_STRESS_RDQ20_MF_H + +#include "active_stress.h" + +/** + * @brief RDQ20-MF mean-field active stress model. + * + * @todo Full model documentation (description, equations, references, units) + * will be added after the complete implementation (increments 3-5). + */ +class RDQ20MF : public ActiveStress { +public: + /// Model label, used for factory registration and XML selection. + static inline const std::string label = "RDQ20-MF"; + + /// @name State vector layout + /// @{ + + /// Number of regulatory-unit (RU) probability states (entries 0-15). + static constexpr unsigned int n_ru_states = 16; + + /// Number of crossbridge (XB) moment states (entries 16-19). + static constexpr unsigned int n_xb_states = 4; + + /// Total number of state variables. + static constexpr unsigned int n_state_variables = n_ru_states + n_xb_states; + + /** + * @brief Flat index of the RU probability state P(TL, TC, TR, CC). + * + * Each argument is 0 or 1 and denotes the state of, respectively, the left + * tropomyosin unit, the central tropomyosin unit, the right tropomyosin unit + * and the central troponin (calcium unbound/bound). The ordering matches the + * reference implementation's serialization (TL outermost, CC innermost) and + * spans [0, 15]. + */ + static constexpr int ru_index(int TL, int TC, int TR, int CC) { + return 8 * TL + 4 * TC + 2 * TR + CC; + } + + /// Flat index of the XB moment state @p i (in [0, 3]), spanning [16, 19]. + static constexpr int xb_index(int i) { + return static_cast(n_ru_states) + i; + } + + /// @} + + /** + * @brief Model parameters class. + * + * @todo Model-specific parameters are added together with the RU/XB dynamics + * in later increments. + */ + class Parameters : public ActiveStressModelParameters { + public: + Parameters() : ActiveStressModelParameters(label) {} + }; + + /** + * @brief Constructor. + */ + RDQ20MF() : ActiveStress(n_state_variables) {} + + /** + * @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. + * + * Sets the state to (1, 0, ..., 0), i.e. all probability mass in the RU state + * P(0, 0, 0, 0) and all crossbridge moments equal to zero. + * + * @param[out] state State vector for a single node, to be initialized by + * this function. + */ + virtual void init_local(Vector &state) const override; + + /** + * @brief Advance in time for a single node. + * + * @note Not implemented yet: the RU and XB state updates are added in later + * increments. Calling this method throws an exception. + */ + 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 active tension for a single node. + * + * @note Temporary: returns zero until the active-tension formula is + * implemented in a later increment. + */ + virtual double + compute_active_tension_local(const Vector &state, + const double fiber_stretch) const override; +}; + +#endif From 64d621f06fadab075e48c70a6bc17512c78b7666 Mon Sep 17 00:00:00 2001 From: Sami Bismar Date: Thu, 23 Jul 2026 19:38:54 -0500 Subject: [PATCH 03/17] Implement RDQ20-MF regulatory unit dynamics --- Code/Source/solver/active_stress_rdq20_mf.cpp | 162 +++++++++++++++++- Code/Source/solver/active_stress_rdq20_mf.h | 90 +++++++++- 2 files changed, 240 insertions(+), 12 deletions(-) diff --git a/Code/Source/solver/active_stress_rdq20_mf.cpp b/Code/Source/solver/active_stress_rdq20_mf.cpp index c5cbb9155..4be2cba81 100644 --- a/Code/Source/solver/active_stress_rdq20_mf.cpp +++ b/Code/Source/solver/active_stress_rdq20_mf.cpp @@ -3,18 +3,31 @@ #include "active_stress_rdq20_mf.h" -#include "Core/Exception.h" +#include +#include void RDQ20MF::read_model_specific_parameters( const ActiveStressModelParameters ¶ms) { - // TODO: read model-specific parameters (added with the RU/XB dynamics in - // later increments). No parameters to read yet. + Kbasic = params.get_scalar("Kbasic"); + Koff = params.get_scalar("Koff"); + Q = params.get_scalar("Q"); + mu = params.get_scalar("mu"); + gamma = params.get_scalar("gamma"); + Kd0 = params.get_scalar("Kd0"); + alphaKd = params.get_scalar("alphaKd"); + SL0 = params.get_scalar("SL0"); } void RDQ20MF::distribute_model_specific_parameters(const CmMod &cm_mod, const cmType &cm) { - // TODO: distribute model-specific parameters (added with the RU/XB dynamics - // in later increments). No parameters to distribute yet. + cm.bcast(cm_mod, &Kbasic); + cm.bcast(cm_mod, &Koff); + cm.bcast(cm_mod, &Q); + cm.bcast(cm_mod, &mu); + cm.bcast(cm_mod, &gamma); + cm.bcast(cm_mod, &Kd0); + cm.bcast(cm_mod, &alphaKd); + cm.bcast(cm_mod, &SL0); } void RDQ20MF::init_local(Vector &state) const { @@ -29,10 +42,52 @@ void RDQ20MF::advance_time_step_local(const double t, const double dt, const double fiber_stretch, const double fiber_stretch_rate, Vector &state) const { - svmp::not_implemented( - "RDQ20-MF active stress dynamics are not implemented yet: the " - "regulatory-unit and crossbridge state updates are added in a later " - "increment."); + // Convert the svMultiPhysics electromechanics inputs to the reference units. + const double calcium_microM = calcium * calcium_mM_to_microM; + const double dt_seconds = dt * time_ms_to_s; + const double sarcomere_length = SL0 * fiber_stretch; + + // Calcium/stretch-independent central-tropomyosin transition rates. + double rates_T[2][2][2][2]; + ru_transition_rates_tropomyosin(rates_T); + + // Troponin transition rates rates_C[CC][TC]: the calcium-binding row (CC = 0) + // depends on calcium and sarcomere length; the unbinding row (CC = 1) does + // not. + const double calcium_on_rate = + Koff / + (Kd0 - alphaKd * (kd_reference_sarcomere_length - sarcomere_length)) * + calcium_microM; + double rates_C[2][2]; + rates_C[0][0] = calcium_on_rate; + rates_C[0][1] = calcium_on_rate; + rates_C[1][0] = Koff; + rates_C[1][1] = Koff / mu; + + // Deserialize the 16 RU probabilities (entries 0-15). The crossbridge moments + // (entries 16-19) are left untouched by this increment. + double state_RU[2][2][2][2]; + for (int TL = 0; TL < 2; ++TL) + for (int TC = 0; TC < 2; ++TC) + for (int TR = 0; TR < 2; ++TR) + for (int CC = 0; CC < 2; ++CC) + state_RU[TL][TC][TR][CC] = state[ru_index(TL, TC, TR, CC)]; + + // Forward-Euler substepping over the outer time step. The final substep is + // shortened so that the outer step is covered exactly. + double time_advanced = 0.0; + while (time_advanced <= dt_seconds - 1.0e-10) { + const double substep = std::min(ru_substep, dt_seconds - time_advanced); + ru_forward_euler_substep(substep, rates_T, rates_C, state_RU); + time_advanced += substep; + } + + // Serialize the updated RU probabilities back into the state vector. + for (int TL = 0; TL < 2; ++TL) + for (int TC = 0; TC < 2; ++TC) + for (int TR = 0; TR < 2; ++TR) + for (int CC = 0; CC < 2; ++CC) + state[ru_index(TL, TC, TR, CC)] = state_RU[TL][TC][TR][CC]; } double RDQ20MF::compute_active_tension_local(const Vector &state, @@ -42,4 +97,93 @@ double RDQ20MF::compute_active_tension_local(const Vector &state, return 0.0; } +void RDQ20MF::ru_transition_rates_tropomyosin( + double (&rates_T)[2][2][2][2]) const { + for (int TL = 0; TL < 2; ++TL) + for (int TR = 0; TR < 2; ++TR) { + const int permissive_neighbors = TL + TR; + + // Rate of leaving the permissive central state (TC = 1). + const double closing_rate = + Kbasic * std::pow(gamma, 2 - permissive_neighbors); + // Rate of leaving the non-permissive central state (TC = 0). + const double opening_rate = + Q * Kbasic * std::pow(gamma, permissive_neighbors); + + rates_T[TL][1][TR][0] = closing_rate; + rates_T[TL][1][TR][1] = closing_rate; + rates_T[TL][0][TR][0] = opening_rate / mu; + rates_T[TL][0][TR][1] = opening_rate; + } +} + +void RDQ20MF::ru_forward_euler_substep( + double dt, const double (&rates_T)[2][2][2][2], + const double (&rates_C)[2][2], double (&state_RU)[2][2][2][2]) const { + // Probability fluxes from central-unit transitions. + double flux_TC[2][2][2][2]; // central tropomyosin + double flux_CC[2][2][2][2]; // central troponin + for (int TL = 0; TL < 2; ++TL) + for (int TC = 0; TC < 2; ++TC) + for (int TR = 0; TR < 2; ++TR) + for (int CC = 0; CC < 2; ++CC) { + flux_TC[TL][TC][TR][CC] = + state_RU[TL][TC][TR][CC] * rates_T[TL][TC][TR][CC]; + flux_CC[TL][TC][TR][CC] = + state_RU[TL][TC][TR][CC] * rates_C[CC][TC]; + } + + // Effective transition rates of the boundary neighbours, obtained from the + // mean-field closure by conditioning the central-unit flux on the neighbour + // pair state. + double rate_left[2][2]; + double rate_right[2][2]; + for (int TL = 0; TL < 2; ++TL) + for (int TC = 0; TC < 2; ++TC) { + double flux_sum = 0.0; + double prob_sum = 0.0; + for (int TR = 0; TR < 2; ++TR) + for (int CC = 0; CC < 2; ++CC) { + flux_sum += flux_TC[TL][TC][TR][CC]; + prob_sum += state_RU[TL][TC][TR][CC]; + } + rate_left[TL][TC] = (prob_sum > 1.0e-12) ? flux_sum / prob_sum : 0.0; + } + for (int TR = 0; TR < 2; ++TR) + for (int TC = 0; TC < 2; ++TC) { + double flux_sum = 0.0; + double prob_sum = 0.0; + for (int TL = 0; TL < 2; ++TL) + for (int CC = 0; CC < 2; ++CC) { + flux_sum += flux_TC[TL][TC][TR][CC]; + prob_sum += state_RU[TL][TC][TR][CC]; + } + rate_right[TR][TC] = (prob_sum > 1.0e-12) ? flux_sum / prob_sum : 0.0; + } + + // Probability fluxes from the boundary-neighbour transitions. + double flux_TL[2][2][2][2]; // left tropomyosin + double flux_TR[2][2][2][2]; // right tropomyosin + for (int TL = 0; TL < 2; ++TL) + for (int TC = 0; TC < 2; ++TC) + for (int TR = 0; TR < 2; ++TR) + for (int CC = 0; CC < 2; ++CC) { + flux_TL[TL][TC][TR][CC] = + state_RU[TL][TC][TR][CC] * rate_left[TC][TL]; + flux_TR[TL][TC][TR][CC] = + state_RU[TL][TC][TR][CC] * rate_right[TC][TR]; + } + + // Forward-Euler update of the 16 RU probabilities. + for (int TL = 0; TL < 2; ++TL) + for (int TC = 0; TC < 2; ++TC) + for (int TR = 0; TR < 2; ++TR) + for (int CC = 0; CC < 2; ++CC) + state_RU[TL][TC][TR][CC] += + dt * (-flux_TL[TL][TC][TR][CC] + flux_TL[1 - TL][TC][TR][CC] - + flux_TC[TL][TC][TR][CC] + flux_TC[TL][1 - TC][TR][CC] - + flux_TR[TL][TC][TR][CC] + flux_TR[TL][TC][1 - TR][CC] - + flux_CC[TL][TC][TR][CC] + flux_CC[TL][TC][TR][1 - CC]); +} + REGISTER_ACTIVE_STRESS_MODEL("RDQ20-MF", RDQ20MF); diff --git a/Code/Source/solver/active_stress_rdq20_mf.h b/Code/Source/solver/active_stress_rdq20_mf.h index e4f557baf..680aa7172 100644 --- a/Code/Source/solver/active_stress_rdq20_mf.h +++ b/Code/Source/solver/active_stress_rdq20_mf.h @@ -52,12 +52,25 @@ class RDQ20MF : public ActiveStress { /** * @brief Model parameters class. * - * @todo Model-specific parameters are added together with the RU/XB dynamics - * in later increments. + * Declares the regulatory-unit (RU) parameters. Crossbridge and active-tension + * parameters are added in later increments. The default values are the human + * body-temperature calibration of the reference implementation, expressed in + * the reference's own units (see the corresponding members below). */ class Parameters : public ActiveStressModelParameters { public: - Parameters() : ActiveStressModelParameters(label) {} + Parameters() : ActiveStressModelParameters(label) { + constexpr bool required = true; + + add_parameter("Kbasic", 13.0, required); + add_parameter("Koff", 100.0, required); + add_parameter("Q", 2.0, required); + add_parameter("mu", 10.0, required); + add_parameter("gamma", 12.0, required); + add_parameter("Kd0", 0.381, required); + add_parameter("alphaKd", -0.571, required); + add_parameter("SL0", 2.2, required); + } }; /** @@ -118,6 +131,77 @@ class RDQ20MF : public ActiveStress { virtual double compute_active_tension_local(const Vector &state, const double fiber_stretch) const override; + +private: + /// @name Regulatory-unit (RU) dynamics helpers + /// @{ + + /** + * @brief Compute the central-tropomyosin transition rate for each local RU + * configuration. + * + * Fills @p rates_T, indexed as @c rates_T[TL][TC][TR][CC], where TL, TC and TR + * are the binary left, central and right tropomyosin states and CC is the + * central troponin calcium-binding state. Each entry is the rate at which the + * central tropomyosin changes state for that configuration. Because the rate + * depends on the neighbour states TL and TR, nearest-neighbour cooperativity + * is retained through the tracked TL-TC-TR configuration. These rates depend + * only on the model parameters, not on calcium or stretch. + */ + void ru_transition_rates_tropomyosin(double (&rates_T)[2][2][2][2]) const; + + /** + * @brief Advance the 16 RU-state probabilities by one forward-Euler substep. + * + * Computes the probability fluxes caused by central-state transitions and the + * effective boundary-neighbour transitions from the mean-field closure, then + * updates @p state_RU in place. + * + * @param[in] dt Substep size [s]. + * @param[in] rates_T Central-tropomyosin transition rates, + * indexed @c rates_T[TL][TC][TR][CC]. + * @param[in] rates_C Troponin transition rates, indexed @c rates_C[CC][TC]. + * @param[in,out] state_RU The 16 RU-state probabilities, + * indexed @c state_RU[TL][TC][TR][CC]. + */ + void ru_forward_euler_substep(double dt, + const double (&rates_T)[2][2][2][2], + const double (&rates_C)[2][2], + double (&state_RU)[2][2][2][2]) const; + + /// @} + + /// @name Interface unit conversions (svMultiPhysics EM units to reference units) + /// @{ + + /// Calcium conversion, millimolar [mM] to micromolar [uM]. + static constexpr double calcium_mM_to_microM = 1.0e3; + + /// Time conversion, milliseconds [ms] to seconds [s]. + static constexpr double time_ms_to_s = 1.0e-3; + + /// @} + + /// RU forward-Euler substep size [s]. + static constexpr double ru_substep = 2.5e-5; + + /// Fixed reference sarcomere length [um] in the length-dependent dissociation + /// constant (distinct from the parameter SL0). + static constexpr double kd_reference_sarcomere_length = 2.15; + + /// @name RU model parameters (reference units) + /// @{ + + double Kbasic; ///< Basic tropomyosin transition rate [s^-1]. + double Koff; ///< Troponin unbinding rate [s^-1]. + double Q; ///< Tropomyosin transition-rate asymmetry factor [-]. + double mu; ///< Calcium-binding cooperativity factor [-]. + double gamma; ///< Nearest-neighbour cooperativity factor [-]. + double Kd0; ///< Calcium dissociation constant at reference length [uM]. + double alphaKd; ///< Length dependence of the dissociation constant [uM/um]. + double SL0; ///< Reference sarcomere length [um]; maps stretch to length. + + /// @} }; #endif From bd776157fa4fab522c7f88e5bb7f5f28676f7b39 Mon Sep 17 00:00:00 2001 From: Sami Bismar Date: Thu, 23 Jul 2026 23:46:07 -0500 Subject: [PATCH 04/17] Implement RDQ20-MF crossbridge dynamics --- Code/Source/solver/active_stress_rdq20_mf.cpp | 84 +++++++++++++++++++ Code/Source/solver/active_stress_rdq20_mf.h | 40 ++++++++- 2 files changed, 120 insertions(+), 4 deletions(-) diff --git a/Code/Source/solver/active_stress_rdq20_mf.cpp b/Code/Source/solver/active_stress_rdq20_mf.cpp index 4be2cba81..7fc1c6a9b 100644 --- a/Code/Source/solver/active_stress_rdq20_mf.cpp +++ b/Code/Source/solver/active_stress_rdq20_mf.cpp @@ -3,6 +3,8 @@ #include "active_stress_rdq20_mf.h" +#include "eigen3/Eigen/Dense" + #include #include @@ -16,6 +18,11 @@ void RDQ20MF::read_model_specific_parameters( Kd0 = params.get_scalar("Kd0"); alphaKd = params.get_scalar("alphaKd"); SL0 = params.get_scalar("SL0"); + + r0 = params.get_scalar("r0"); + alpha = params.get_scalar("alpha"); + mu0_fP = params.get_scalar("mu0_fP"); + mu1_fP = params.get_scalar("mu1_fP"); } void RDQ20MF::distribute_model_specific_parameters(const CmMod &cm_mod, @@ -28,6 +35,11 @@ void RDQ20MF::distribute_model_specific_parameters(const CmMod &cm_mod, cm.bcast(cm_mod, &Kd0); cm.bcast(cm_mod, &alphaKd); cm.bcast(cm_mod, &SL0); + + cm.bcast(cm_mod, &r0); + cm.bcast(cm_mod, &alpha); + cm.bcast(cm_mod, &mu0_fP); + cm.bcast(cm_mod, &mu1_fP); } void RDQ20MF::init_local(Vector &state) const { @@ -82,12 +94,25 @@ void RDQ20MF::advance_time_step_local(const double t, const double dt, time_advanced += substep; } + // Advance the crossbridge moments (entries 16-19) from the updated RU state. + // The reference velocity v = -dSL/dt / SL0 reduces to -d(lambda)/dt because + // SL = SL0 * lambda; here it is expressed in reference time units [s^-1]. + const double velocity = -fiber_stretch_rate / time_ms_to_s; + double state_XB[4]; + for (int i = 0; i < 4; ++i) + state_XB[i] = state[xb_index(i)]; + xb_implicit_update(dt_seconds, velocity, rates_T, state_RU, state_XB); + // Serialize the updated RU probabilities back into the state vector. for (int TL = 0; TL < 2; ++TL) for (int TC = 0; TC < 2; ++TC) for (int TR = 0; TR < 2; ++TR) for (int CC = 0; CC < 2; ++CC) state[ru_index(TL, TC, TR, CC)] = state_RU[TL][TC][TR][CC]; + + // Serialize the updated crossbridge moments back into the state vector. + for (int i = 0; i < 4; ++i) + state[xb_index(i)] = state_XB[i]; } double RDQ20MF::compute_active_tension_local(const Vector &state, @@ -186,4 +211,63 @@ void RDQ20MF::ru_forward_euler_substep( flux_CC[TL][TC][TR][CC] + flux_CC[TL][TC][TR][1 - CC]); } +void RDQ20MF::xb_implicit_update(double dt, double velocity, + const double (&rates_T)[2][2][2][2], + const double (&state_RU)[2][2][2][2], + double (&state_XB)[4]) const { + // Permissivity and the permissive/non-permissive probability fluxes from the + // updated RU state. + double permissivity = 0.0; + double flux_PN = 0.0; + double flux_NP = 0.0; + for (int TL = 0; TL < 2; ++TL) + for (int TR = 0; TR < 2; ++TR) + for (int CC = 0; CC < 2; ++CC) { + permissivity += state_RU[TL][1][TR][CC]; + flux_PN += state_RU[TL][1][TR][CC] * rates_T[TL][1][TR][CC]; + flux_NP += state_RU[TL][0][TR][CC] * rates_T[TL][0][TR][CC]; + } + + // Effective permissive->non-permissive and non-permissive->permissive rates. + const double k_PN = (permissivity >= 1.0e-12) ? flux_PN / permissivity : 0.0; + const double k_NP = + ((1.0 - permissivity) >= 1.0e-12) ? flux_NP / (1.0 - permissivity) : 0.0; + + // Use the calibrated RDQ20-MF specialization of the general XB system: + // new XBs attach only in the permissive state (f_N = 0), and both XB + // populations share r(v) = r0 + alpha * |v|. Non-permissive moments + // are populated by P-to-N transitions of already-attached XBs. + const double r = r0 + alpha * std::abs(velocity); + const double diag_P = r + k_PN; + const double diag_N = r + k_NP; + + // Implicit-Euler system (I - dt * A) x = rhs for the four moments. The matrix + // is zero-initialized so the structurally-zero entries are correct. + Eigen::Matrix system = Eigen::Matrix::Zero(); + system(0, 0) = -diag_P; + system(1, 1) = -diag_P; + system(2, 2) = -diag_N; + system(3, 3) = -diag_N; + system(0, 2) = k_NP; + system(1, 3) = k_NP; + system(2, 0) = k_PN; + system(3, 1) = k_PN; + system(1, 0) = -velocity; + system(3, 2) = -velocity; + system *= -dt; + for (int i = 0; i < 4; ++i) + system(i, i) += 1.0; + + Eigen::Matrix rhs; + rhs(0) = state_XB[0] + dt * permissivity * mu0_fP; + rhs(1) = state_XB[1] + dt * permissivity * mu1_fP; + rhs(2) = state_XB[2]; + rhs(3) = state_XB[3]; + + const Eigen::Matrix solution = + system.colPivHouseholderQr().solve(rhs); + for (int i = 0; i < 4; ++i) + state_XB[i] = solution(i); +} + REGISTER_ACTIVE_STRESS_MODEL("RDQ20-MF", RDQ20MF); diff --git a/Code/Source/solver/active_stress_rdq20_mf.h b/Code/Source/solver/active_stress_rdq20_mf.h index 680aa7172..c76958562 100644 --- a/Code/Source/solver/active_stress_rdq20_mf.h +++ b/Code/Source/solver/active_stress_rdq20_mf.h @@ -52,10 +52,11 @@ class RDQ20MF : public ActiveStress { /** * @brief Model parameters class. * - * Declares the regulatory-unit (RU) parameters. Crossbridge and active-tension - * parameters are added in later increments. The default values are the human - * body-temperature calibration of the reference implementation, expressed in - * the reference's own units (see the corresponding members below). + * Declares the regulatory-unit (RU) and crossbridge (XB) parameters. The + * active-tension parameters are added in a later increment. The default values + * are the human body-temperature calibration of the reference implementation, + * expressed in the reference's own units (see the corresponding members + * below). */ class Parameters : public ActiveStressModelParameters { public: @@ -70,6 +71,11 @@ class RDQ20MF : public ActiveStress { add_parameter("Kd0", 0.381, required); add_parameter("alphaKd", -0.571, required); add_parameter("SL0", 2.2, required); + + add_parameter("r0", 134.31, required); + add_parameter("alpha", 25.184, required); + add_parameter("mu0_fP", 32.653, required); + add_parameter("mu1_fP", 0.778, required); } }; @@ -169,6 +175,27 @@ class RDQ20MF : public ActiveStress { const double (&rates_C)[2][2], double (&state_RU)[2][2][2][2]) const; + /** + * @brief Advance the four crossbridge moments by one implicit-Euler step. + * + * Computes the permissivity and the effective permissive/non-permissive + * transition rates from the updated RU probabilities, forms the 4x4 linear + * system for the implicit update and solves it in place for @p state_XB. + * + * @param[in] dt Outer time step [s]. + * @param[in] velocity Shortening velocity @f$-\dot{SL}/SL_0@f$ [s^-1]. + * @param[in] rates_T Central-tropomyosin transition rates, + * indexed @c rates_T[TL][TC][TR][CC]. + * @param[in] state_RU The updated 16 RU-state probabilities, + * indexed @c state_RU[TL][TC][TR][CC]. + * @param[in,out] state_XB The four crossbridge moments, ordered + * @f$[\mu_P^0, \mu_P^1, \mu_N^0, \mu_N^1]@f$. + */ + void xb_implicit_update(double dt, double velocity, + const double (&rates_T)[2][2][2][2], + const double (&state_RU)[2][2][2][2], + double (&state_XB)[4]) const; + /// @} /// @name Interface unit conversions (svMultiPhysics EM units to reference units) @@ -201,6 +228,11 @@ class RDQ20MF : public ActiveStress { double alphaKd; ///< Length dependence of the dissociation constant [uM/um]. double SL0; ///< Reference sarcomere length [um]; maps stretch to length. + double r0; ///< Combined attachment-detachment rate at zero velocity [s^-1]. + double alpha; ///< Coefficient of |v| in r(v) = r0 + alpha * |v| [-]. + double mu0_fP; ///< Permissive influx into the zeroth-moment crossbridge state [s^-1]. + double mu1_fP; ///< Permissive influx into the first-moment crossbridge state [s^-1]. + /// @} }; From ee3cf55685cb719483fca0a852b47f48cb960dbc Mon Sep 17 00:00:00 2001 From: Sami Bismar Date: Fri, 24 Jul 2026 12:05:08 -0500 Subject: [PATCH 05/17] Implement RDQ20-MF active tension --- Code/Source/solver/active_stress_rdq20_mf.cpp | 36 ++++++- Code/Source/solver/active_stress_rdq20_mf.h | 93 ++++++++++++++++--- 2 files changed, 115 insertions(+), 14 deletions(-) diff --git a/Code/Source/solver/active_stress_rdq20_mf.cpp b/Code/Source/solver/active_stress_rdq20_mf.cpp index 7fc1c6a9b..c66697347 100644 --- a/Code/Source/solver/active_stress_rdq20_mf.cpp +++ b/Code/Source/solver/active_stress_rdq20_mf.cpp @@ -23,6 +23,11 @@ void RDQ20MF::read_model_specific_parameters( alpha = params.get_scalar("alpha"); mu0_fP = params.get_scalar("mu0_fP"); mu1_fP = params.get_scalar("mu1_fP"); + + LA = params.get_scalar("LA"); + LM = params.get_scalar("LM"); + LB = params.get_scalar("LB"); + a_XB = params.get_scalar("a_XB"); } void RDQ20MF::distribute_model_specific_parameters(const CmMod &cm_mod, @@ -40,6 +45,11 @@ void RDQ20MF::distribute_model_specific_parameters(const CmMod &cm_mod, cm.bcast(cm_mod, &alpha); cm.bcast(cm_mod, &mu0_fP); cm.bcast(cm_mod, &mu1_fP); + + cm.bcast(cm_mod, &LA); + cm.bcast(cm_mod, &LM); + cm.bcast(cm_mod, &LB); + cm.bcast(cm_mod, &a_XB); } void RDQ20MF::init_local(Vector &state) const { @@ -117,9 +127,14 @@ void RDQ20MF::advance_time_step_local(const double t, const double dt, double RDQ20MF::compute_active_tension_local(const Vector &state, const double fiber_stretch) const { - // TEMPORARY: returns zero until the active-tension formula is implemented in - // a later increment. - return 0.0; + const double sarcomere_length = SL0 * fiber_stretch; + + // Active tension from the permissive and non-permissive XB first moments + // (state entries 17 and 19), scaled by the single-overlap fraction and the + // upscaling factor a_XB. The moments and overlap fraction are dimensionless, + // so the returned tension has the same stress unit as a_XB (no conversion). + return a_XB * (state[xb_index(1)] + state[xb_index(3)]) * + fraction_single_overlap(sarcomere_length); } void RDQ20MF::ru_transition_rates_tropomyosin( @@ -270,4 +285,19 @@ void RDQ20MF::xb_implicit_update(double dt, double velocity, state_XB[i] = solution(i); } +double RDQ20MF::fraction_single_overlap(double sarcomere_length) const { + const double SL = sarcomere_length; + const double half_single_overlap = (LM - LB) * 0.5; + + if (SL > LA && SL <= LM) + return (SL - LA) / half_single_overlap; + if (SL > LM && SL <= 2.0 * LA - LB) + return (SL + LM - 2.0 * LA) * 0.5 / half_single_overlap; + if (SL > 2.0 * LA - LB && SL <= 2.0 * LA + LB) + return 1.0; + if (SL > 2.0 * LA + LB && SL <= 2.0 * LA + LM) + return (LM + 2.0 * LA - SL) * 0.5 / half_single_overlap; + return 0.0; +} + REGISTER_ACTIVE_STRESS_MODEL("RDQ20-MF", RDQ20MF); diff --git a/Code/Source/solver/active_stress_rdq20_mf.h b/Code/Source/solver/active_stress_rdq20_mf.h index c76958562..d0073442a 100644 --- a/Code/Source/solver/active_stress_rdq20_mf.h +++ b/Code/Source/solver/active_stress_rdq20_mf.h @@ -9,8 +9,36 @@ /** * @brief RDQ20-MF mean-field active stress model. * - * @todo Full model documentation (description, equations, references, units) - * will be added after the complete implementation (increments 3-5). + * This class implements the mean-field RDQ20-MF sarcomere model of cardiomyocyte + * force generation described in [1] and is validated against the authors' + * reference implementation [2]. The node-local state has 20 variables: 16 + * regulatory-unit (RU) probabilities (entries 0-15) describing the + * tropomyosin/troponin configuration of a triplet of neighbouring units, and 4 + * crossbridge (XB) moments (entries 16-19). The RU probabilities are advanced + * with an explicit forward-Euler substepping scheme and the XB moments with one + * implicit-Euler step per time step; the active tension is then reconstructed + * from the XB first moments. + * + * The active tension is + * @f[ + * T_\text{act} = a_\text{XB} \, (\mu_P^1 + \mu_N^1) \, \phi(SL)\;, + * @f] + * where @f$\mu_P^1@f$ and @f$\mu_N^1@f$ are the permissive and non-permissive + * first XB moments (state entries 17 and 19), @f$\phi(SL)@f$ is the single-overlap + * fraction of the sarcomere at sarcomere length @f$SL = SL_0 \, \lambda@f$ (with + * @f$\lambda@f$ the fiber stretch), and @f$a_\text{XB}@f$ is the tension + * upscaling factor. Because @f$\mu_P^1 + \mu_N^1@f$ and @f$\phi(SL)@f$ are + * dimensionless, @f$a_\text{XB}@f$ sets the units of the returned active tension; + * it is stored in the stress units of the simulation, so no separate output + * conversion is applied. + * + * The model is calibrated in a fixed unit system (calcium in [uM], time in [s], + * length in [um]); the svMultiPhysics inputs are converted to these units at the + * interface (see the conversion members below). + * + * **References**: + * 1. [Regazzoni, Dede', Quarteroni (2020)](https://doi.org/10.1371/journal.pcbi.1008294) + * 2. [F. Regazzoni, cardiac-activation reference implementation](https://github.com/FrancescoRegazzoni/cardiac-activation) */ class RDQ20MF : public ActiveStress { public: @@ -52,11 +80,18 @@ class RDQ20MF : public ActiveStress { /** * @brief Model parameters class. * - * Declares the regulatory-unit (RU) and crossbridge (XB) parameters. The - * active-tension parameters are added in a later increment. The default values - * are the human body-temperature calibration of the reference implementation, - * expressed in the reference's own units (see the corresponding members - * below). + * Declares the regulatory-unit (RU), crossbridge (XB), geometry, and tension + * parameters. The values registered below correspond to the published human + * body-temperature calibration of the reference implementation, expressed in + * the units documented for each member. The registered value of @c a_XB, + * 22.894, expresses the reference calibration in MPa; the value supplied in + * solver.xml must instead use the stress unit of the simulation's mechanical + * configuration. + * + * All parameters are required. A complete @c RDQ20-MF parameter block containing + * every parameter must be provided in solver.xml. Any value may be changed to + * use a different calibration, but omitting a parameter causes a parse error; + * the registered reference value is not used as an automatic default. */ class Parameters : public ActiveStressModelParameters { public: @@ -76,6 +111,11 @@ class RDQ20MF : public ActiveStress { add_parameter("alpha", 25.184, required); add_parameter("mu0_fP", 32.653, required); add_parameter("mu1_fP", 0.778, required); + + add_parameter("LA", 1.25, required); + add_parameter("LM", 1.65, required); + add_parameter("LB", 0.18, required); + add_parameter("a_XB", 22.894, required); } }; @@ -119,8 +159,10 @@ class RDQ20MF : public ActiveStress { /** * @brief Advance in time for a single node. * - * @note Not implemented yet: the RU and XB state updates are added in later - * increments. Calling this method throws an exception. + * Advances the RU probabilities (entries 0-15) with the forward-Euler + * substepping scheme and then the XB moments (entries 16-19) with one + * implicit-Euler step, using the calcium, fiber stretch and fiber-stretch rate + * at the node. */ virtual void advance_time_step_local(const double t, const double dt, const double calcium, @@ -131,8 +173,10 @@ class RDQ20MF : public ActiveStress { /** * @brief Compute the active tension for a single node. * - * @note Temporary: returns zero until the active-tension formula is - * implemented in a later increment. + * Evaluates @f$T_\text{act} = a_\text{XB} (\mu_P^1 + \mu_N^1) \phi(SL)@f$ from + * the XB first moments (state entries 17 and 19) and the single-overlap + * fraction at @f$SL = SL_0 \, \lambda@f$, where @f$\lambda@f$ is @p + * fiber_stretch. The result is in the stress units of the simulation. */ virtual double compute_active_tension_local(const Vector &state, @@ -196,6 +240,17 @@ class RDQ20MF : public ActiveStress { const double (&state_RU)[2][2][2][2], double (&state_XB)[4]) const; + /** + * @brief Single-overlap fraction of the sarcomere at a given length. + * + * Returns the fraction @f$\phi(SL) \in [0, 1]@f$ of the sarcomere over which + * thin and thick filaments overlap exactly once, a piecewise-linear function + * of the sarcomere length built from the filament geometry (LA, LM, LB). + * + * @param[in] sarcomere_length Sarcomere length @f$SL@f$ [um]. + */ + double fraction_single_overlap(double sarcomere_length) const; + /// @} /// @name Interface unit conversions (svMultiPhysics EM units to reference units) @@ -233,6 +288,22 @@ class RDQ20MF : public ActiveStress { double mu0_fP; ///< Permissive influx into the zeroth-moment crossbridge state [s^-1]. double mu1_fP; ///< Permissive influx into the first-moment crossbridge state [s^-1]. + double LA; ///< Thin-filament (actin) length [um]. + double LM; ///< Thick-filament (myosin) length [um]. + double LB; ///< Length of the myosin bare zone [um]. + + /// Tension upscaling factor. + /// + /// Because the crossbridge moments and the overlap fraction are dimensionless, + /// a_XB is the only quantity carrying stress units, so the returned active + /// tension has the same stress unit as a_XB and no stress-unit conversion is + /// performed. a_XB must therefore be expressed in the same stress unit as the + /// mechanical configuration. The reference calibration value 22.894 is + /// expressed in MPa, consistent with the coupled electromechanics slab case; + /// the equivalent values are 22.894e3 in kPa and 22.894e6 in Pa. Provide the + /// value matching the case's stress unit. + double a_XB; + /// @} }; From acebfec5626a8a16384bb1409730114a59bcfb82 Mon Sep 17 00:00:00 2001 From: Sami Bismar Date: Sat, 25 Jul 2026 12:40:14 -0500 Subject: [PATCH 06/17] Add RDQ20-MF electromechanics regression test --- .../electromechanics/slab_RDQ20MF/README.md | 84 ++++++++ .../electromechanics/slab_RDQ20MF/mesh/X0.vtp | 3 + .../electromechanics/slab_RDQ20MF/mesh/X1.vtp | 3 + .../slab_RDQ20MF/mesh/volume.vtu | 3 + .../slab_RDQ20MF/result_001.vtu | 3 + .../electromechanics/slab_RDQ20MF/solver.xml | 195 ++++++++++++++++++ tests/test_electromechanics.py | 5 + 7 files changed, 296 insertions(+) create mode 100755 tests/cases/electromechanics/slab_RDQ20MF/README.md create mode 100644 tests/cases/electromechanics/slab_RDQ20MF/mesh/X0.vtp create mode 100644 tests/cases/electromechanics/slab_RDQ20MF/mesh/X1.vtp create mode 100644 tests/cases/electromechanics/slab_RDQ20MF/mesh/volume.vtu create mode 100644 tests/cases/electromechanics/slab_RDQ20MF/result_001.vtu create mode 100644 tests/cases/electromechanics/slab_RDQ20MF/solver.xml diff --git a/tests/cases/electromechanics/slab_RDQ20MF/README.md b/tests/cases/electromechanics/slab_RDQ20MF/README.md new file mode 100755 index 000000000..44dae8ea6 --- /dev/null +++ b/tests/cases/electromechanics/slab_RDQ20MF/README.md @@ -0,0 +1,84 @@ + +# **Problem Description** + +Simulate cardiac electromechanics on a slab of myocardial tissue. This test +couples 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. + +## Electrophysiology + +The propagation of the transmembrane potential is modeled with the +ten-Tusscher-Panfilov (`TTP`) cell activation model [2, 3], using epicardial +parameters (included from `../../cep/ttp_parameters/ttp_epicardium_parameters.xml`) +and anisotropic conductivity aligned with the fiber direction. The domain is split +into two `Domain`s: an unstimulated region (`domain 1`) and a stimulated region +(`domain 2`) where an external `Istim` stimulus initiates depolarization. + +``` + + -35.714 + 0.0 + 2.0 + 10000.0 + +``` + +## Mechanics + +The tissue is modeled as a nearly incompressible Holzapfel-Ogden material with +modified anisotropy (`HolzapfelOgden-ModifiedAnisotropy`) [4]. Active contraction +is driven by the calcium concentration computed by the electrophysiology model, +through the RDQ20-MF mean-field active-stress model [5], configured with the +published human body-temperature calibration. The scalar active tension is +distributed along the fiber, sheet, and sheet-normal directions using the same +directional weights as the NashPanfilov slab case. + +``` + + RDQ20-MF + + 0.7 + 0.2 + 0.1 + + ... + +``` + +The slab is fixed with a zero-displacement Dirichlet boundary condition on the +`X1` face, and contracts as the depolarization wave propagates through the tissue. + +## Regression Reference + +The RDQ20-MF scalar active tension in `result_001.vtu` -- reconstructed as +`Active_tension_fibers + Active_tension_sheets + Active_tension_normal` (the +directional weights sum to one) -- was compared node-by-node against the authors' +RDQ20-MF reference implementation at commit +`26f05df28891df7b3c69f16bb136cdced6b63c4d`, using the calcium and stretch inputs +from this one-step test. The results agreed to machine precision. The remaining +fields in the VTU serve as integrated svMultiPhysics regression references and +were not independently validated by the RDQ20-MF reference code. + +## References + +[1] S. A. Niederer, E. Kerfoot, A. P. Benson, et al. Verification of cardiac tissue +electrophysiology simulators using an N-version benchmark. Philosophical Transactions +of the Royal Society A, 369(1954):4331–4351, 2011. + +[2] K. H. W. J. ten Tusscher, D. Noble, P. J. Noble, and A. V. Panfilov. A model for +human ventricular tissue. American Journal of Physiology-Heart and Circulatory +Physiology, 286(4):H1573–H1589, apr 2004. + +[3] K. H. W. J. ten Tusscher and A. V. Panfilov. Alternans and spiral breakup in a +human ventricular tissue model. American Journal of Physiology-Heart and Circulatory +Physiology, 291(3):H1088–H1100, sep 2006. + +[4] G. A. Holzapfel and R. W. Ogden. Constitutive modelling of passive myocardium: a +structurally based framework for material characterization. Philosophical Transactions +of the Royal Society A, 367(1902):3445–3475, 2009. + +[5] F. Regazzoni, L. Dede', and A. Quarteroni. Biophysically detailed mathematical +models of multiscale cardiac active mechanics. PLOS Computational Biology, +16(10):e1008294, 2020. diff --git a/tests/cases/electromechanics/slab_RDQ20MF/mesh/X0.vtp b/tests/cases/electromechanics/slab_RDQ20MF/mesh/X0.vtp new file mode 100644 index 000000000..eaebcc0eb --- /dev/null +++ b/tests/cases/electromechanics/slab_RDQ20MF/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_RDQ20MF/mesh/X1.vtp b/tests/cases/electromechanics/slab_RDQ20MF/mesh/X1.vtp new file mode 100644 index 000000000..07ad82435 --- /dev/null +++ b/tests/cases/electromechanics/slab_RDQ20MF/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_RDQ20MF/mesh/volume.vtu b/tests/cases/electromechanics/slab_RDQ20MF/mesh/volume.vtu new file mode 100644 index 000000000..dd12aca89 --- /dev/null +++ b/tests/cases/electromechanics/slab_RDQ20MF/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_RDQ20MF/result_001.vtu b/tests/cases/electromechanics/slab_RDQ20MF/result_001.vtu new file mode 100644 index 000000000..fc338238d --- /dev/null +++ b/tests/cases/electromechanics/slab_RDQ20MF/result_001.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:31c6a24f61d33d7646970715f4e5c0dc359cd9b2fad077bdcfbd7625f428bbec +size 1443224 diff --git a/tests/cases/electromechanics/slab_RDQ20MF/solver.xml b/tests/cases/electromechanics/slab_RDQ20MF/solver.xml new file mode 100644 index 000000000..8a0071b60 --- /dev/null +++ b/tests/cases/electromechanics/slab_RDQ20MF/solver.xml @@ -0,0 +1,195 @@ + + + + false + 3 + 1 + 1.0 + 0.50 + STOP_SIM + + true + result + 1 + 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-12 + + + 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-12 + 50 + + + + + 1 + 6 + 1e-12 + + 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 + + + + RDQ20-MF + + + 0.7 + 0.2 + 0.1 + + + + 13.0 + 100.0 + 2.0 + 10.0 + 12.0 + 0.381 + -0.571 + 2.2 + 134.31 + 25.184 + 32.653 + 0.778 + 1.25 + 1.65 + 0.18 + 22.894 + + + + + true + true + true + true + true + true + true + true + + true + true + true + + + + + fsils + + 1e-12 + 1e-14 + 1000 + + + + Dir + 0.0 + + + + + diff --git a/tests/test_electromechanics.py b/tests/test_electromechanics.py index 5eaa0cdad..b27f5e48f 100644 --- a/tests/test_electromechanics.py +++ b/tests/test_electromechanics.py @@ -26,3 +26,8 @@ def test_slab(n_proc): test_folder = "slab" run_with_reference(base_folder, test_folder, fields, n_proc, t_max=1) + + +def test_slab_RDQ20MF(n_proc): + test_folder = "slab_RDQ20MF" + run_with_reference(base_folder, test_folder, fields, n_proc, t_max=1) From 1b2830066c59335716cbb153e4d6115f4b7981fa Mon Sep 17 00:00:00 2001 From: Sami Bismar Date: Tue, 28 Jul 2026 10:17:18 -0500 Subject: [PATCH 07/17] Fix RDQ20MF rate naming, alphaKd guard, and minor cleanups --- Code/Source/solver/active_stress_rdq20_mf.cpp | 16 +++++++++++++--- .../electromechanics/slab_RDQ20MF/README.md | 0 2 files changed, 13 insertions(+), 3 deletions(-) mode change 100755 => 100644 tests/cases/electromechanics/slab_RDQ20MF/README.md diff --git a/Code/Source/solver/active_stress_rdq20_mf.cpp b/Code/Source/solver/active_stress_rdq20_mf.cpp index c66697347..54e9f13e1 100644 --- a/Code/Source/solver/active_stress_rdq20_mf.cpp +++ b/Code/Source/solver/active_stress_rdq20_mf.cpp @@ -17,6 +17,12 @@ void RDQ20MF::read_model_specific_parameters( gamma = params.get_scalar("gamma"); Kd0 = params.get_scalar("Kd0"); alphaKd = params.get_scalar("alphaKd"); + if (alphaKd > 0.0) + svmp::raise( + "RDQ20MF: alphaKd must be <= 0 (positive values reduce calcium " + "sensitivity with stretch, reversing length-dependent activation, " + "and can produce a zero dissociation constant at physiological " + "sarcomere lengths)."); SL0 = params.get_scalar("SL0"); r0 = params.get_scalar("r0"); @@ -53,7 +59,7 @@ void RDQ20MF::distribute_model_specific_parameters(const CmMod &cm_mod, } void RDQ20MF::init_local(Vector &state) const { - for (unsigned int i = 0; i < n_states; ++i) + for (unsigned int i = 0; i < n_state_variables; ++i) state[i] = 0.0; state[ru_index(0, 0, 0, 0)] = 1.0; // == state[0] @@ -202,6 +208,10 @@ void RDQ20MF::ru_forward_euler_substep( } // Probability fluxes from the boundary-neighbour transitions. + // TL's only neighbour is TC on its right → rate_right[TC][TL]. + // TR's only neighbour is TC on its left → rate_left[TC][TR]. + // (rate_left == rate_right numerically due to mean-field LR symmetry, so the + // result is unchanged, but the names now match the physical convention.) double flux_TL[2][2][2][2]; // left tropomyosin double flux_TR[2][2][2][2]; // right tropomyosin for (int TL = 0; TL < 2; ++TL) @@ -209,9 +219,9 @@ void RDQ20MF::ru_forward_euler_substep( for (int TR = 0; TR < 2; ++TR) for (int CC = 0; CC < 2; ++CC) { flux_TL[TL][TC][TR][CC] = - state_RU[TL][TC][TR][CC] * rate_left[TC][TL]; + state_RU[TL][TC][TR][CC] * rate_right[TC][TL]; flux_TR[TL][TC][TR][CC] = - state_RU[TL][TC][TR][CC] * rate_right[TC][TR]; + state_RU[TL][TC][TR][CC] * rate_left[TC][TR]; } // Forward-Euler update of the 16 RU probabilities. diff --git a/tests/cases/electromechanics/slab_RDQ20MF/README.md b/tests/cases/electromechanics/slab_RDQ20MF/README.md old mode 100755 new mode 100644 From 920769676fff1fff904238e35671b6e191eb3721 Mon Sep 17 00:00:00 2001 From: Sami Bismar Date: Fri, 31 Jul 2026 12:16:14 -0500 Subject: [PATCH 08/17] Normalize RDQ20-MF active tension for PK2 assembly --- Code/Source/solver/active_stress_rdq20_mf.cpp | 18 ++++++++--- Code/Source/solver/active_stress_rdq20_mf.h | 27 +++++++++++++--- .../electromechanics/slab_RDQ20MF/README.md | 31 ++++++++++++++----- 3 files changed, 58 insertions(+), 18 deletions(-) diff --git a/Code/Source/solver/active_stress_rdq20_mf.cpp b/Code/Source/solver/active_stress_rdq20_mf.cpp index 54e9f13e1..3ab9db87b 100644 --- a/Code/Source/solver/active_stress_rdq20_mf.cpp +++ b/Code/Source/solver/active_stress_rdq20_mf.cpp @@ -3,7 +3,9 @@ #include "active_stress_rdq20_mf.h" +#include "FE/Common/FEException.h" #include "eigen3/Eigen/Dense" +#include "utils.h" #include #include @@ -133,14 +135,20 @@ void RDQ20MF::advance_time_step_local(const double t, const double dt, double RDQ20MF::compute_active_tension_local(const Vector &state, const double fiber_stretch) const { + if (utils::is_zero(fiber_stretch)) + svmp::raise( + "RDQ20MF: fiber_stretch is zero or near zero; this indicates a " + "degenerate (collapsed) element and is not a valid deformation state."); + const double sarcomere_length = SL0 * fiber_stretch; - // Active tension from the permissive and non-permissive XB first moments - // (state entries 17 and 19), scaled by the single-overlap fraction and the - // upscaling factor a_XB. The moments and overlap fraction are dimensionless, - // so the returned tension has the same stress unit as a_XB (no conversion). + // Compute paper tension T̃_act = a_XB * (μ_P^1 + μ_N^1) * φ(SL) from the + // permissive and non-permissive XB first moments (state entries 17 and 19), + // scaled by the single-overlap fraction and the upscaling factor a_XB. + // Return T̃_act / λ_f so that the assembled first-PK active stress has norm + // T̃_act, consistent with the RDQ20-MF paper formulation. return a_XB * (state[xb_index(1)] + state[xb_index(3)]) * - fraction_single_overlap(sarcomere_length); + fraction_single_overlap(sarcomere_length) / fiber_stretch; } void RDQ20MF::ru_transition_rates_tropomyosin( diff --git a/Code/Source/solver/active_stress_rdq20_mf.h b/Code/Source/solver/active_stress_rdq20_mf.h index d0073442a..57e3ec6db 100644 --- a/Code/Source/solver/active_stress_rdq20_mf.h +++ b/Code/Source/solver/active_stress_rdq20_mf.h @@ -171,12 +171,29 @@ class RDQ20MF : public ActiveStress { Vector &state) const override; /** - * @brief Compute the active tension for a single node. + * @brief Compute the PK2 active-tension coefficient for a single node. * - * Evaluates @f$T_\text{act} = a_\text{XB} (\mu_P^1 + \mu_N^1) \phi(SL)@f$ from - * the XB first moments (state entries 17 and 19) and the single-overlap - * fraction at @f$SL = SL_0 \, \lambda@f$, where @f$\lambda@f$ is @p - * fiber_stretch. The result is in the stress units of the simulation. + * Computes the paper tension + * @f$\widetilde{T}_\text{act} = a_\text{XB} (\mu_P^1 + \mu_N^1) \phi(SL)@f$ + * from the XB first moments (state entries 17 and 19) and the single-overlap + * fraction at @f$SL = SL_0 \lambda_f@f$, then returns + * @f$\widetilde{T}_\text{act} / \lambda_f@f$. + * + * The division by @f$\lambda_f@f$ is necessary because svMultiPhysics + * assembles the active stress as + * @f$\mathbf{S}_\text{act} = T_\text{code}\,\mathbf{f}_0\otimes\mathbf{f}_0@f$ + * (second-PK), so the resulting first-PK active stress satisfies + * @f$\|\mathbf{P}_\text{act}\| = \widetilde{T}_\text{act}@f$, matching the + * paper's definition of active tension. + * + * The returned value has the stress units of @f$a_\text{XB}@f$. + * + * @note The slab regression test uses a (0.7, 0.2, 0.1) fiber/sheet/normal + * directional distribution, which is an svMultiPhysics extension of the + * paper's fiber-only active stress formulation. + * + * @throws svmp::FE::InvalidArgumentException if @p fiber_stretch is zero + * or near zero (degenerate/collapsed element). */ virtual double compute_active_tension_local(const Vector &state, diff --git a/tests/cases/electromechanics/slab_RDQ20MF/README.md b/tests/cases/electromechanics/slab_RDQ20MF/README.md index 44dae8ea6..70dcc6b18 100644 --- a/tests/cases/electromechanics/slab_RDQ20MF/README.md +++ b/tests/cases/electromechanics/slab_RDQ20MF/README.md @@ -50,16 +50,31 @@ directional weights as the NashPanfilov slab case. The slab is fixed with a zero-displacement Dirichlet boundary condition on the `X1` face, and contracts as the depolarization wave propagates through the tissue. +The (0.7, 0.2, 0.1) fiber/sheet/sheet-normal directional weights are an +svMultiPhysics extension of the paper's fiber-only active stress formulation; they +are not prescribed by the RDQ20-MF model itself. + ## Regression Reference -The RDQ20-MF scalar active tension in `result_001.vtu` -- reconstructed as -`Active_tension_fibers + Active_tension_sheets + Active_tension_normal` (the -directional weights sum to one) -- was compared node-by-node against the authors' -RDQ20-MF reference implementation at commit -`26f05df28891df7b3c69f16bb136cdced6b63c4d`, using the calcium and stretch inputs -from this one-step test. The results agreed to machine precision. The remaining -fields in the VTU serve as integrated svMultiPhysics regression references and -were not independently validated by the RDQ20-MF reference code. +svMultiPhysics stores `T_code = T̃_act / λ_f` in each per-direction field +(`Active_tension_fibers`, `Active_tension_sheets`, `Active_tension_normal`), where +`T̃_act` is the paper tension from Regazzoni et al. (2020) and `λ_f` is the fiber +stretch; the paper tension is recovered as `T̃_act = λ_f · T_code`. Each +per-direction field stores `η · T_code`, where `η` is the corresponding directional +weight; their sum recovers `T_code`. + +The active tension fields in `result_001.vtu` were validated node-by-node against +the C++ reference implementation at commit +[`26f05df`](https://github.com/FrancescoRegazzoni/cardiac-activation/commit/26f05df28891df7b3c69f16bb136cdced6b63c4d). +Both implementations use the same implicit-Euler XB scheme, so agreement is +to machine precision (~1e-16 relative error). The comparison evaluates +`λ_f · T_code` (i.e. `T̃_act`) from the svMultiPhysics output against the reference +C++ active tension, using the calcium and sarcomere-length inputs from this one-step +test. Because this test covers only one time step from the undeformed configuration, +`λ_f ≈ 1.0` throughout, and the regression reference does not strongly distinguish +the normalized (`T̃_act / λ_f`) and unnormalized (`T̃_act`) conventions. +The remaining fields in the VTU serve as integrated svMultiPhysics regression +references and were not independently validated by the RDQ20-MF reference code. ## References From 7f34f444a198eeceec9a48e101c04474d64268e4 Mon Sep 17 00:00:00 2001 From: Sami Bismar Date: Fri, 31 Jul 2026 20:18:11 -0500 Subject: [PATCH 09/17] Rename RDQ20MF model to MeanFieldActiveStress --- Code/Source/solver/CMakeLists.txt | 2 +- ...20_mf.cpp => active_stress_mean_field.cpp} | 44 +++++++++---------- ..._rdq20_mf.h => active_stress_mean_field.h} | 20 ++++----- .../README.md | 2 +- .../mesh/X0.vtp | 0 .../mesh/X1.vtp | 0 .../mesh/volume.vtu | 0 .../result_001.vtu | 0 .../solver.xml | 6 +-- tests/test_electromechanics.py | 4 +- 10 files changed, 39 insertions(+), 39 deletions(-) rename Code/Source/solver/{active_stress_rdq20_mf.cpp => active_stress_mean_field.cpp} (88%) rename Code/Source/solver/{active_stress_rdq20_mf.h => active_stress_mean_field.h} (94%) rename tests/cases/electromechanics/{slab_RDQ20MF => slab_MeanFieldActiveStress}/README.md (99%) rename tests/cases/electromechanics/{slab_RDQ20MF => slab_MeanFieldActiveStress}/mesh/X0.vtp (100%) rename tests/cases/electromechanics/{slab_RDQ20MF => slab_MeanFieldActiveStress}/mesh/X1.vtp (100%) rename tests/cases/electromechanics/{slab_RDQ20MF => slab_MeanFieldActiveStress}/mesh/volume.vtu (100%) rename tests/cases/electromechanics/{slab_RDQ20MF => slab_MeanFieldActiveStress}/result_001.vtu (100%) rename tests/cases/electromechanics/{slab_RDQ20MF => slab_MeanFieldActiveStress}/solver.xml (98%) diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt index f4688f25d..f4f04d347 100644 --- a/Code/Source/solver/CMakeLists.txt +++ b/Code/Source/solver/CMakeLists.txt @@ -234,7 +234,7 @@ set(CSRCS active_stress_uniform_unsteady.cpp active_stress_ode.cpp active_stress_nash_panfilov.cpp - active_stress_rdq20_mf.cpp + active_stress_mean_field.cpp SPLIT.c diff --git a/Code/Source/solver/active_stress_rdq20_mf.cpp b/Code/Source/solver/active_stress_mean_field.cpp similarity index 88% rename from Code/Source/solver/active_stress_rdq20_mf.cpp rename to Code/Source/solver/active_stress_mean_field.cpp index 3ab9db87b..48f6ac5d4 100644 --- a/Code/Source/solver/active_stress_rdq20_mf.cpp +++ b/Code/Source/solver/active_stress_mean_field.cpp @@ -1,7 +1,7 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the // University of California, and others. SPDX-License-Identifier: BSD-3-Clause -#include "active_stress_rdq20_mf.h" +#include "active_stress_mean_field.h" #include "FE/Common/FEException.h" #include "eigen3/Eigen/Dense" @@ -10,7 +10,7 @@ #include #include -void RDQ20MF::read_model_specific_parameters( +void MeanFieldActiveStress::read_model_specific_parameters( const ActiveStressModelParameters ¶ms) { Kbasic = params.get_scalar("Kbasic"); Koff = params.get_scalar("Koff"); @@ -21,7 +21,7 @@ void RDQ20MF::read_model_specific_parameters( alphaKd = params.get_scalar("alphaKd"); if (alphaKd > 0.0) svmp::raise( - "RDQ20MF: alphaKd must be <= 0 (positive values reduce calcium " + "MeanFieldActiveStress: alphaKd must be <= 0 (positive values reduce calcium " "sensitivity with stretch, reversing length-dependent activation, " "and can produce a zero dissociation constant at physiological " "sarcomere lengths)."); @@ -38,8 +38,8 @@ void RDQ20MF::read_model_specific_parameters( a_XB = params.get_scalar("a_XB"); } -void RDQ20MF::distribute_model_specific_parameters(const CmMod &cm_mod, - const cmType &cm) { +void MeanFieldActiveStress::distribute_model_specific_parameters( + const CmMod &cm_mod, const cmType &cm) { cm.bcast(cm_mod, &Kbasic); cm.bcast(cm_mod, &Koff); cm.bcast(cm_mod, &Q); @@ -60,18 +60,17 @@ void RDQ20MF::distribute_model_specific_parameters(const CmMod &cm_mod, cm.bcast(cm_mod, &a_XB); } -void RDQ20MF::init_local(Vector &state) const { +void MeanFieldActiveStress::init_local(Vector &state) const { for (unsigned int i = 0; i < n_state_variables; ++i) state[i] = 0.0; state[ru_index(0, 0, 0, 0)] = 1.0; // == state[0] } -void RDQ20MF::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 { +void MeanFieldActiveStress::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 { // Convert the svMultiPhysics electromechanics inputs to the reference units. const double calcium_microM = calcium * calcium_mM_to_microM; const double dt_seconds = dt * time_ms_to_s; @@ -133,11 +132,11 @@ void RDQ20MF::advance_time_step_local(const double t, const double dt, state[xb_index(i)] = state_XB[i]; } -double RDQ20MF::compute_active_tension_local(const Vector &state, - const double fiber_stretch) const { +double MeanFieldActiveStress::compute_active_tension_local( + const Vector &state, const double fiber_stretch) const { if (utils::is_zero(fiber_stretch)) svmp::raise( - "RDQ20MF: fiber_stretch is zero or near zero; this indicates a " + "MeanFieldActiveStress: fiber_stretch is zero or near zero; this indicates a " "degenerate (collapsed) element and is not a valid deformation state."); const double sarcomere_length = SL0 * fiber_stretch; @@ -151,7 +150,7 @@ double RDQ20MF::compute_active_tension_local(const Vector &state, fraction_single_overlap(sarcomere_length) / fiber_stretch; } -void RDQ20MF::ru_transition_rates_tropomyosin( +void MeanFieldActiveStress::ru_transition_rates_tropomyosin( double (&rates_T)[2][2][2][2]) const { for (int TL = 0; TL < 2; ++TL) for (int TR = 0; TR < 2; ++TR) { @@ -171,7 +170,7 @@ void RDQ20MF::ru_transition_rates_tropomyosin( } } -void RDQ20MF::ru_forward_euler_substep( +void MeanFieldActiveStress::ru_forward_euler_substep( double dt, const double (&rates_T)[2][2][2][2], const double (&rates_C)[2][2], double (&state_RU)[2][2][2][2]) const { // Probability fluxes from central-unit transitions. @@ -244,10 +243,11 @@ void RDQ20MF::ru_forward_euler_substep( flux_CC[TL][TC][TR][CC] + flux_CC[TL][TC][TR][1 - CC]); } -void RDQ20MF::xb_implicit_update(double dt, double velocity, - const double (&rates_T)[2][2][2][2], - const double (&state_RU)[2][2][2][2], - double (&state_XB)[4]) const { +void MeanFieldActiveStress::xb_implicit_update( + double dt, double velocity, + const double (&rates_T)[2][2][2][2], + const double (&state_RU)[2][2][2][2], + double (&state_XB)[4]) const { // Permissivity and the permissive/non-permissive probability fluxes from the // updated RU state. double permissivity = 0.0; @@ -303,7 +303,7 @@ void RDQ20MF::xb_implicit_update(double dt, double velocity, state_XB[i] = solution(i); } -double RDQ20MF::fraction_single_overlap(double sarcomere_length) const { +double MeanFieldActiveStress::fraction_single_overlap(double sarcomere_length) const { const double SL = sarcomere_length; const double half_single_overlap = (LM - LB) * 0.5; @@ -318,4 +318,4 @@ double RDQ20MF::fraction_single_overlap(double sarcomere_length) const { return 0.0; } -REGISTER_ACTIVE_STRESS_MODEL("RDQ20-MF", RDQ20MF); +REGISTER_ACTIVE_STRESS_MODEL("MeanFieldActiveStress", MeanFieldActiveStress); diff --git a/Code/Source/solver/active_stress_rdq20_mf.h b/Code/Source/solver/active_stress_mean_field.h similarity index 94% rename from Code/Source/solver/active_stress_rdq20_mf.h rename to Code/Source/solver/active_stress_mean_field.h index 57e3ec6db..29b5a1294 100644 --- a/Code/Source/solver/active_stress_rdq20_mf.h +++ b/Code/Source/solver/active_stress_mean_field.h @@ -1,17 +1,17 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the // University of California, and others. SPDX-License-Identifier: BSD-3-Clause -#ifndef ACTIVE_STRESS_RDQ20_MF_H -#define ACTIVE_STRESS_RDQ20_MF_H +#ifndef ACTIVE_STRESS_MEAN_FIELD_H +#define ACTIVE_STRESS_MEAN_FIELD_H #include "active_stress.h" /** - * @brief RDQ20-MF mean-field active stress model. + * @brief Mean-field active stress model (implements the RDQ20-MF formulation). * * This class implements the mean-field RDQ20-MF sarcomere model of cardiomyocyte - * force generation described in [1] and is validated against the authors' - * reference implementation [2]. The node-local state has 20 variables: 16 + * force generation of Regazzoni, Dedè, and Quarteroni (2020), described in [1] + * and validated against the authors' reference implementation [2]. The node-local state has 20 variables: 16 * regulatory-unit (RU) probabilities (entries 0-15) describing the * tropomyosin/troponin configuration of a triplet of neighbouring units, and 4 * crossbridge (XB) moments (entries 16-19). The RU probabilities are advanced @@ -40,10 +40,10 @@ * 1. [Regazzoni, Dede', Quarteroni (2020)](https://doi.org/10.1371/journal.pcbi.1008294) * 2. [F. Regazzoni, cardiac-activation reference implementation](https://github.com/FrancescoRegazzoni/cardiac-activation) */ -class RDQ20MF : public ActiveStress { +class MeanFieldActiveStress : public ActiveStress { public: /// Model label, used for factory registration and XML selection. - static inline const std::string label = "RDQ20-MF"; + static inline const std::string label = "MeanFieldActiveStress"; /// @name State vector layout /// @{ @@ -88,8 +88,8 @@ class RDQ20MF : public ActiveStress { * solver.xml must instead use the stress unit of the simulation's mechanical * configuration. * - * All parameters are required. A complete @c RDQ20-MF parameter block containing - * every parameter must be provided in solver.xml. Any value may be changed to + * All parameters are required. A complete @c MeanFieldActiveStress parameter block + * containing every parameter must be provided in solver.xml. Any value may be changed to * use a different calibration, but omitting a parameter causes a parse error; * the registered reference value is not used as an automatic default. */ @@ -122,7 +122,7 @@ class RDQ20MF : public ActiveStress { /** * @brief Constructor. */ - RDQ20MF() : ActiveStress(n_state_variables) {} + MeanFieldActiveStress() : ActiveStress(n_state_variables) {} /** * @brief Construct an instance of model parameters. diff --git a/tests/cases/electromechanics/slab_RDQ20MF/README.md b/tests/cases/electromechanics/slab_MeanFieldActiveStress/README.md similarity index 99% rename from tests/cases/electromechanics/slab_RDQ20MF/README.md rename to tests/cases/electromechanics/slab_MeanFieldActiveStress/README.md index 70dcc6b18..d416b7396 100644 --- a/tests/cases/electromechanics/slab_RDQ20MF/README.md +++ b/tests/cases/electromechanics/slab_MeanFieldActiveStress/README.md @@ -37,7 +37,7 @@ directional weights as the NashPanfilov slab case. ``` - RDQ20-MF + MeanFieldActiveStress 0.7 0.2 diff --git a/tests/cases/electromechanics/slab_RDQ20MF/mesh/X0.vtp b/tests/cases/electromechanics/slab_MeanFieldActiveStress/mesh/X0.vtp similarity index 100% rename from tests/cases/electromechanics/slab_RDQ20MF/mesh/X0.vtp rename to tests/cases/electromechanics/slab_MeanFieldActiveStress/mesh/X0.vtp diff --git a/tests/cases/electromechanics/slab_RDQ20MF/mesh/X1.vtp b/tests/cases/electromechanics/slab_MeanFieldActiveStress/mesh/X1.vtp similarity index 100% rename from tests/cases/electromechanics/slab_RDQ20MF/mesh/X1.vtp rename to tests/cases/electromechanics/slab_MeanFieldActiveStress/mesh/X1.vtp diff --git a/tests/cases/electromechanics/slab_RDQ20MF/mesh/volume.vtu b/tests/cases/electromechanics/slab_MeanFieldActiveStress/mesh/volume.vtu similarity index 100% rename from tests/cases/electromechanics/slab_RDQ20MF/mesh/volume.vtu rename to tests/cases/electromechanics/slab_MeanFieldActiveStress/mesh/volume.vtu diff --git a/tests/cases/electromechanics/slab_RDQ20MF/result_001.vtu b/tests/cases/electromechanics/slab_MeanFieldActiveStress/result_001.vtu similarity index 100% rename from tests/cases/electromechanics/slab_RDQ20MF/result_001.vtu rename to tests/cases/electromechanics/slab_MeanFieldActiveStress/result_001.vtu diff --git a/tests/cases/electromechanics/slab_RDQ20MF/solver.xml b/tests/cases/electromechanics/slab_MeanFieldActiveStress/solver.xml similarity index 98% rename from tests/cases/electromechanics/slab_RDQ20MF/solver.xml rename to tests/cases/electromechanics/slab_MeanFieldActiveStress/solver.xml index 8a0071b60..2f67d9cc5 100644 --- a/tests/cases/electromechanics/slab_RDQ20MF/solver.xml +++ b/tests/cases/electromechanics/slab_MeanFieldActiveStress/solver.xml @@ -133,7 +133,7 @@ - RDQ20-MF + MeanFieldActiveStress 0.7 @@ -141,7 +141,7 @@ 0.1 - + 13.0 100.0 2.0 @@ -158,7 +158,7 @@ 1.65 0.18 22.894 - + diff --git a/tests/test_electromechanics.py b/tests/test_electromechanics.py index b27f5e48f..ed3b7b642 100644 --- a/tests/test_electromechanics.py +++ b/tests/test_electromechanics.py @@ -28,6 +28,6 @@ def test_slab(n_proc): run_with_reference(base_folder, test_folder, fields, n_proc, t_max=1) -def test_slab_RDQ20MF(n_proc): - test_folder = "slab_RDQ20MF" +def test_slab_MeanFieldActiveStress(n_proc): + test_folder = "slab_MeanFieldActiveStress" run_with_reference(base_folder, test_folder, fields, n_proc, t_max=1) From a59075497902ac7565ca42514d500efa484d7a93 Mon Sep 17 00:00:00 2001 From: Sami Bismar Date: Tue, 4 Aug 2026 12:40:59 -0500 Subject: [PATCH 10/17] Rename mean-field active stress model to Regazzoni --- Code/Source/solver/CMakeLists.txt | 2 +- ..._field.cpp => active_stress_regazzoni.cpp} | 26 +++++++++---------- ...mean_field.h => active_stress_regazzoni.h} | 12 ++++----- .../README.md | 2 +- .../mesh/X0.vtp | 0 .../mesh/X1.vtp | 0 .../mesh/volume.vtu | 0 .../result_001.vtu | 0 .../solver.xml | 6 ++--- tests/test_electromechanics.py | 4 +-- 10 files changed, 26 insertions(+), 26 deletions(-) rename Code/Source/solver/{active_stress_mean_field.cpp => active_stress_regazzoni.cpp} (93%) rename Code/Source/solver/{active_stress_mean_field.h => active_stress_regazzoni.h} (97%) rename tests/cases/electromechanics/{slab_MeanFieldActiveStress => slab_Regazzoni}/README.md (99%) rename tests/cases/electromechanics/{slab_MeanFieldActiveStress => slab_Regazzoni}/mesh/X0.vtp (100%) rename tests/cases/electromechanics/{slab_MeanFieldActiveStress => slab_Regazzoni}/mesh/X1.vtp (100%) rename tests/cases/electromechanics/{slab_MeanFieldActiveStress => slab_Regazzoni}/mesh/volume.vtu (100%) rename tests/cases/electromechanics/{slab_MeanFieldActiveStress => slab_Regazzoni}/result_001.vtu (100%) rename tests/cases/electromechanics/{slab_MeanFieldActiveStress => slab_Regazzoni}/solver.xml (98%) diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt index f4f04d347..b37e8b985 100644 --- a/Code/Source/solver/CMakeLists.txt +++ b/Code/Source/solver/CMakeLists.txt @@ -234,7 +234,7 @@ set(CSRCS active_stress_uniform_unsteady.cpp active_stress_ode.cpp active_stress_nash_panfilov.cpp - active_stress_mean_field.cpp + active_stress_regazzoni.cpp SPLIT.c diff --git a/Code/Source/solver/active_stress_mean_field.cpp b/Code/Source/solver/active_stress_regazzoni.cpp similarity index 93% rename from Code/Source/solver/active_stress_mean_field.cpp rename to Code/Source/solver/active_stress_regazzoni.cpp index 48f6ac5d4..32a4290d1 100644 --- a/Code/Source/solver/active_stress_mean_field.cpp +++ b/Code/Source/solver/active_stress_regazzoni.cpp @@ -1,7 +1,7 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the // University of California, and others. SPDX-License-Identifier: BSD-3-Clause -#include "active_stress_mean_field.h" +#include "active_stress_regazzoni.h" #include "FE/Common/FEException.h" #include "eigen3/Eigen/Dense" @@ -10,7 +10,7 @@ #include #include -void MeanFieldActiveStress::read_model_specific_parameters( +void RegazzoniActiveStress::read_model_specific_parameters( const ActiveStressModelParameters ¶ms) { Kbasic = params.get_scalar("Kbasic"); Koff = params.get_scalar("Koff"); @@ -21,7 +21,7 @@ void MeanFieldActiveStress::read_model_specific_parameters( alphaKd = params.get_scalar("alphaKd"); if (alphaKd > 0.0) svmp::raise( - "MeanFieldActiveStress: alphaKd must be <= 0 (positive values reduce calcium " + "RegazzoniActiveStress: alphaKd must be <= 0 (positive values reduce calcium " "sensitivity with stretch, reversing length-dependent activation, " "and can produce a zero dissociation constant at physiological " "sarcomere lengths)."); @@ -38,7 +38,7 @@ void MeanFieldActiveStress::read_model_specific_parameters( a_XB = params.get_scalar("a_XB"); } -void MeanFieldActiveStress::distribute_model_specific_parameters( +void RegazzoniActiveStress::distribute_model_specific_parameters( const CmMod &cm_mod, const cmType &cm) { cm.bcast(cm_mod, &Kbasic); cm.bcast(cm_mod, &Koff); @@ -60,14 +60,14 @@ void MeanFieldActiveStress::distribute_model_specific_parameters( cm.bcast(cm_mod, &a_XB); } -void MeanFieldActiveStress::init_local(Vector &state) const { +void RegazzoniActiveStress::init_local(Vector &state) const { for (unsigned int i = 0; i < n_state_variables; ++i) state[i] = 0.0; state[ru_index(0, 0, 0, 0)] = 1.0; // == state[0] } -void MeanFieldActiveStress::advance_time_step_local( +void RegazzoniActiveStress::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 { @@ -132,11 +132,11 @@ void MeanFieldActiveStress::advance_time_step_local( state[xb_index(i)] = state_XB[i]; } -double MeanFieldActiveStress::compute_active_tension_local( +double RegazzoniActiveStress::compute_active_tension_local( const Vector &state, const double fiber_stretch) const { if (utils::is_zero(fiber_stretch)) svmp::raise( - "MeanFieldActiveStress: fiber_stretch is zero or near zero; this indicates a " + "RegazzoniActiveStress: fiber_stretch is zero or near zero; this indicates a " "degenerate (collapsed) element and is not a valid deformation state."); const double sarcomere_length = SL0 * fiber_stretch; @@ -150,7 +150,7 @@ double MeanFieldActiveStress::compute_active_tension_local( fraction_single_overlap(sarcomere_length) / fiber_stretch; } -void MeanFieldActiveStress::ru_transition_rates_tropomyosin( +void RegazzoniActiveStress::ru_transition_rates_tropomyosin( double (&rates_T)[2][2][2][2]) const { for (int TL = 0; TL < 2; ++TL) for (int TR = 0; TR < 2; ++TR) { @@ -170,7 +170,7 @@ void MeanFieldActiveStress::ru_transition_rates_tropomyosin( } } -void MeanFieldActiveStress::ru_forward_euler_substep( +void RegazzoniActiveStress::ru_forward_euler_substep( double dt, const double (&rates_T)[2][2][2][2], const double (&rates_C)[2][2], double (&state_RU)[2][2][2][2]) const { // Probability fluxes from central-unit transitions. @@ -243,7 +243,7 @@ void MeanFieldActiveStress::ru_forward_euler_substep( flux_CC[TL][TC][TR][CC] + flux_CC[TL][TC][TR][1 - CC]); } -void MeanFieldActiveStress::xb_implicit_update( +void RegazzoniActiveStress::xb_implicit_update( double dt, double velocity, const double (&rates_T)[2][2][2][2], const double (&state_RU)[2][2][2][2], @@ -303,7 +303,7 @@ void MeanFieldActiveStress::xb_implicit_update( state_XB[i] = solution(i); } -double MeanFieldActiveStress::fraction_single_overlap(double sarcomere_length) const { +double RegazzoniActiveStress::fraction_single_overlap(double sarcomere_length) const { const double SL = sarcomere_length; const double half_single_overlap = (LM - LB) * 0.5; @@ -318,4 +318,4 @@ double MeanFieldActiveStress::fraction_single_overlap(double sarcomere_length) c return 0.0; } -REGISTER_ACTIVE_STRESS_MODEL("MeanFieldActiveStress", MeanFieldActiveStress); +REGISTER_ACTIVE_STRESS_MODEL("Regazzoni", RegazzoniActiveStress); diff --git a/Code/Source/solver/active_stress_mean_field.h b/Code/Source/solver/active_stress_regazzoni.h similarity index 97% rename from Code/Source/solver/active_stress_mean_field.h rename to Code/Source/solver/active_stress_regazzoni.h index 29b5a1294..1f67d1aa1 100644 --- a/Code/Source/solver/active_stress_mean_field.h +++ b/Code/Source/solver/active_stress_regazzoni.h @@ -1,8 +1,8 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the // University of California, and others. SPDX-License-Identifier: BSD-3-Clause -#ifndef ACTIVE_STRESS_MEAN_FIELD_H -#define ACTIVE_STRESS_MEAN_FIELD_H +#ifndef ACTIVE_STRESS_REGAZZONI_H +#define ACTIVE_STRESS_REGAZZONI_H #include "active_stress.h" @@ -40,10 +40,10 @@ * 1. [Regazzoni, Dede', Quarteroni (2020)](https://doi.org/10.1371/journal.pcbi.1008294) * 2. [F. Regazzoni, cardiac-activation reference implementation](https://github.com/FrancescoRegazzoni/cardiac-activation) */ -class MeanFieldActiveStress : public ActiveStress { +class RegazzoniActiveStress : public ActiveStress { public: /// Model label, used for factory registration and XML selection. - static inline const std::string label = "MeanFieldActiveStress"; + static inline const std::string label = "Regazzoni"; /// @name State vector layout /// @{ @@ -88,7 +88,7 @@ class MeanFieldActiveStress : public ActiveStress { * solver.xml must instead use the stress unit of the simulation's mechanical * configuration. * - * All parameters are required. A complete @c MeanFieldActiveStress parameter block + * All parameters are required. A complete @c Regazzoni parameter block * containing every parameter must be provided in solver.xml. Any value may be changed to * use a different calibration, but omitting a parameter causes a parse error; * the registered reference value is not used as an automatic default. @@ -122,7 +122,7 @@ class MeanFieldActiveStress : public ActiveStress { /** * @brief Constructor. */ - MeanFieldActiveStress() : ActiveStress(n_state_variables) {} + RegazzoniActiveStress() : ActiveStress(n_state_variables) {} /** * @brief Construct an instance of model parameters. diff --git a/tests/cases/electromechanics/slab_MeanFieldActiveStress/README.md b/tests/cases/electromechanics/slab_Regazzoni/README.md similarity index 99% rename from tests/cases/electromechanics/slab_MeanFieldActiveStress/README.md rename to tests/cases/electromechanics/slab_Regazzoni/README.md index d416b7396..57d8e92c4 100644 --- a/tests/cases/electromechanics/slab_MeanFieldActiveStress/README.md +++ b/tests/cases/electromechanics/slab_Regazzoni/README.md @@ -37,7 +37,7 @@ directional weights as the NashPanfilov slab case. ``` - MeanFieldActiveStress + Regazzoni 0.7 0.2 diff --git a/tests/cases/electromechanics/slab_MeanFieldActiveStress/mesh/X0.vtp b/tests/cases/electromechanics/slab_Regazzoni/mesh/X0.vtp similarity index 100% rename from tests/cases/electromechanics/slab_MeanFieldActiveStress/mesh/X0.vtp rename to tests/cases/electromechanics/slab_Regazzoni/mesh/X0.vtp diff --git a/tests/cases/electromechanics/slab_MeanFieldActiveStress/mesh/X1.vtp b/tests/cases/electromechanics/slab_Regazzoni/mesh/X1.vtp similarity index 100% rename from tests/cases/electromechanics/slab_MeanFieldActiveStress/mesh/X1.vtp rename to tests/cases/electromechanics/slab_Regazzoni/mesh/X1.vtp diff --git a/tests/cases/electromechanics/slab_MeanFieldActiveStress/mesh/volume.vtu b/tests/cases/electromechanics/slab_Regazzoni/mesh/volume.vtu similarity index 100% rename from tests/cases/electromechanics/slab_MeanFieldActiveStress/mesh/volume.vtu rename to tests/cases/electromechanics/slab_Regazzoni/mesh/volume.vtu diff --git a/tests/cases/electromechanics/slab_MeanFieldActiveStress/result_001.vtu b/tests/cases/electromechanics/slab_Regazzoni/result_001.vtu similarity index 100% rename from tests/cases/electromechanics/slab_MeanFieldActiveStress/result_001.vtu rename to tests/cases/electromechanics/slab_Regazzoni/result_001.vtu diff --git a/tests/cases/electromechanics/slab_MeanFieldActiveStress/solver.xml b/tests/cases/electromechanics/slab_Regazzoni/solver.xml similarity index 98% rename from tests/cases/electromechanics/slab_MeanFieldActiveStress/solver.xml rename to tests/cases/electromechanics/slab_Regazzoni/solver.xml index 2f67d9cc5..ddc60007a 100644 --- a/tests/cases/electromechanics/slab_MeanFieldActiveStress/solver.xml +++ b/tests/cases/electromechanics/slab_Regazzoni/solver.xml @@ -133,7 +133,7 @@ - MeanFieldActiveStress + Regazzoni 0.7 @@ -141,7 +141,7 @@ 0.1 - + 13.0 100.0 2.0 @@ -158,7 +158,7 @@ 1.65 0.18 22.894 - + diff --git a/tests/test_electromechanics.py b/tests/test_electromechanics.py index ed3b7b642..92a5d8bb0 100644 --- a/tests/test_electromechanics.py +++ b/tests/test_electromechanics.py @@ -28,6 +28,6 @@ def test_slab(n_proc): run_with_reference(base_folder, test_folder, fields, n_proc, t_max=1) -def test_slab_MeanFieldActiveStress(n_proc): - test_folder = "slab_MeanFieldActiveStress" +def test_slab_regazzoni(n_proc): + test_folder = "slab_Regazzoni" run_with_reference(base_folder, test_folder, fields, n_proc, t_max=1) From ce7b31b72247a8b3a3bcf16f2d473039ed7cd72d Mon Sep 17 00:00:00 2001 From: Sami Bismar Date: Tue, 4 Aug 2026 14:49:40 -0500 Subject: [PATCH 11/17] Update active stress documentation and notation --- Code/Source/solver/active_stress_ode.h | 4 ++-- Code/Source/solver/active_stress_regazzoni.h | 9 +++++---- utilities/fiber_generation/DOCUMENTATION.md | 2 +- 3 files changed, 8 insertions(+), 7 deletions(-) diff --git a/Code/Source/solver/active_stress_ode.h b/Code/Source/solver/active_stress_ode.h index f0696eabe..6ebfa7d09 100644 --- a/Code/Source/solver/active_stress_ode.h +++ b/Code/Source/solver/active_stress_ode.h @@ -16,7 +16,7 @@ * \dv{\astressstate}{t} &= * \mathbf{F}_\text{AS}(t, \astressstate, \calcium, \fiberstretch, * \fiberstretchrate)\;, \\ - * \Tact &= \Tact(\astressstate)\;. + * \Tact &= \Tact(\astressstate, \fiberstretch)\;. * \end{aligned} @f] * * ### Numerical methods @@ -25,7 +25,7 @@ * @ref ODESolver. After that, the active tension is computed for every node * @f$i@f$ as: * @f[ - * {\Tact}_{i}^{n+1} = \Tact(\astressstate_i^{n+1})\;. + * {\Tact}_{i}^{n+1} = \Tact(\astressstate_i^{n+1}, \fiberstretch_i^{n+1})\;. * @f] * * ### Implementing derived models diff --git a/Code/Source/solver/active_stress_regazzoni.h b/Code/Source/solver/active_stress_regazzoni.h index 1f67d1aa1..e152e193e 100644 --- a/Code/Source/solver/active_stress_regazzoni.h +++ b/Code/Source/solver/active_stress_regazzoni.h @@ -10,14 +10,15 @@ * @brief Mean-field active stress model (implements the RDQ20-MF formulation). * * This class implements the mean-field RDQ20-MF sarcomere model of cardiomyocyte - * force generation of Regazzoni, Dedè, and Quarteroni (2020), described in [1] + * force generation of Regazzoni, Dede', and Quarteroni (2020), described in [1] * and validated against the authors' reference implementation [2]. The node-local state has 20 variables: 16 * regulatory-unit (RU) probabilities (entries 0-15) describing the * tropomyosin/troponin configuration of a triplet of neighbouring units, and 4 * crossbridge (XB) moments (entries 16-19). The RU probabilities are advanced - * with an explicit forward-Euler substepping scheme and the XB moments with one - * implicit-Euler step per time step; the active tension is then reconstructed - * from the XB first moments. + * with an explicit forward-Euler substepping scheme — for every macro time step, + * a number of smaller sub-steps are taken to update the RU states — and the XB + * moments with one implicit-Euler step per time step; the active tension is then + * reconstructed from the XB first moments. * * The active tension is * @f[ diff --git a/utilities/fiber_generation/DOCUMENTATION.md b/utilities/fiber_generation/DOCUMENTATION.md index 550c0b129..c9a2dd1c9 100644 --- a/utilities/fiber_generation/DOCUMENTATION.md +++ b/utilities/fiber_generation/DOCUMENTATION.md @@ -539,4 +539,4 @@ For coherency, for all methods and for all chambers, we consider the transmural # References 1. Bayer, J. D., Blake, R. C., Plank, G., & Trayanova, N. A. (2012). A Novel Rule-Based Algorithm for Assigning Myocardial Fiber Orientation to Computational Heart Models. Annals of Biomedical Engineering, 40(10), 2243–2254. https://doi.org/10.1007/s10439-012-0593-5 2. Doste, R., Soto‐Iglesias, D., Bernardino, G., Alcaine, A., Sebastian, R., Giffard‐Roisin, S., Sermesant, M., Berruezo, A., Sanchez‐Quintana, D., & Camara, O. (2019). A rule‐based method to model myocardial fiber orientation in cardiac biventricular geometries with outflow tracts. International Journal for Numerical Methods in Biomedical Engineering, 35(4). https://doi.org/10.1002/cnm.3185 -3. Piersanti, R., Africa, P. C., Fedele, M., Vergara, C., Dedè, L., Corno, A. F., & Quarteroni, A. (2021). Modeling cardiac muscle fibers in ventricular and atrial electrophysiology simulations. Computer Methods in Applied Mechanics and Engineering, 373, 113468. https://doi.org/10.1016/j.cma.2020.113468 \ No newline at end of file +3. Piersanti, R., Africa, P. C., Fedele, M., Vergara, C., Dede', L., Corno, A. F., & Quarteroni, A. (2021). Modeling cardiac muscle fibers in ventricular and atrial electrophysiology simulations. Computer Methods in Applied Mechanics and Engineering, 373, 113468. https://doi.org/10.1016/j.cma.2020.113468 \ No newline at end of file From f0c628ee1ac94cf27edbabae256edc3bfa9d4fd2 Mon Sep 17 00:00:00 2001 From: Sami Bismar Date: Tue, 4 Aug 2026 15:18:30 -0500 Subject: [PATCH 12/17] Use unsigned indices for active stress state helpers --- Code/Source/solver/active_stress_regazzoni.h | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/Code/Source/solver/active_stress_regazzoni.h b/Code/Source/solver/active_stress_regazzoni.h index e152e193e..bfec88529 100644 --- a/Code/Source/solver/active_stress_regazzoni.h +++ b/Code/Source/solver/active_stress_regazzoni.h @@ -67,13 +67,14 @@ class RegazzoniActiveStress : public ActiveStress { * reference implementation's serialization (TL outermost, CC innermost) and * spans [0, 15]. */ - static constexpr int ru_index(int TL, int TC, int TR, int CC) { + static constexpr unsigned int ru_index(unsigned int TL, unsigned int TC, + unsigned int TR, unsigned int CC) { return 8 * TL + 4 * TC + 2 * TR + CC; } /// Flat index of the XB moment state @p i (in [0, 3]), spanning [16, 19]. - static constexpr int xb_index(int i) { - return static_cast(n_ru_states) + i; + static constexpr unsigned int xb_index(unsigned int i) { + return n_ru_states + i; } /// @} From ec1eb58989124167f44cb0509a77e228f4bcd35d Mon Sep 17 00:00:00 2001 From: Sami Bismar Date: Tue, 4 Aug 2026 18:19:56 -0500 Subject: [PATCH 13/17] Remove fixed unit assumptions from Regazzoni active stress model --- .../Source/solver/active_stress_regazzoni.cpp | 16 ++- Code/Source/solver/active_stress_regazzoni.h | 98 +++++++------------ .../slab_Regazzoni/result_001.vtu | 4 +- .../slab_Regazzoni/solver.xml | 32 +++--- 4 files changed, 58 insertions(+), 92 deletions(-) diff --git a/Code/Source/solver/active_stress_regazzoni.cpp b/Code/Source/solver/active_stress_regazzoni.cpp index 32a4290d1..3d8f2a78f 100644 --- a/Code/Source/solver/active_stress_regazzoni.cpp +++ b/Code/Source/solver/active_stress_regazzoni.cpp @@ -71,9 +71,6 @@ void RegazzoniActiveStress::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 { - // Convert the svMultiPhysics electromechanics inputs to the reference units. - const double calcium_microM = calcium * calcium_mM_to_microM; - const double dt_seconds = dt * time_ms_to_s; const double sarcomere_length = SL0 * fiber_stretch; // Calcium/stretch-independent central-tropomyosin transition rates. @@ -86,7 +83,7 @@ void RegazzoniActiveStress::advance_time_step_local( const double calcium_on_rate = Koff / (Kd0 - alphaKd * (kd_reference_sarcomere_length - sarcomere_length)) * - calcium_microM; + calcium; double rates_C[2][2]; rates_C[0][0] = calcium_on_rate; rates_C[0][1] = calcium_on_rate; @@ -105,20 +102,19 @@ void RegazzoniActiveStress::advance_time_step_local( // Forward-Euler substepping over the outer time step. The final substep is // shortened so that the outer step is covered exactly. double time_advanced = 0.0; - while (time_advanced <= dt_seconds - 1.0e-10) { - const double substep = std::min(ru_substep, dt_seconds - time_advanced); + while (time_advanced <= dt - 1.0e-10) { + const double substep = std::min(ru_substep, dt - time_advanced); ru_forward_euler_substep(substep, rates_T, rates_C, state_RU); time_advanced += substep; } // Advance the crossbridge moments (entries 16-19) from the updated RU state. - // The reference velocity v = -dSL/dt / SL0 reduces to -d(lambda)/dt because - // SL = SL0 * lambda; here it is expressed in reference time units [s^-1]. - const double velocity = -fiber_stretch_rate / time_ms_to_s; + // The velocity v = -dSL/dt / SL0 reduces to -d(lambda)/dt because SL = SL0 * lambda. + const double velocity = -fiber_stretch_rate; double state_XB[4]; for (int i = 0; i < 4; ++i) state_XB[i] = state[xb_index(i)]; - xb_implicit_update(dt_seconds, velocity, rates_T, state_RU, state_XB); + xb_implicit_update(dt, velocity, rates_T, state_RU, state_XB); // Serialize the updated RU probabilities back into the state vector. for (int TL = 0; TL < 2; ++TL) diff --git a/Code/Source/solver/active_stress_regazzoni.h b/Code/Source/solver/active_stress_regazzoni.h index bfec88529..be82e8c56 100644 --- a/Code/Source/solver/active_stress_regazzoni.h +++ b/Code/Source/solver/active_stress_regazzoni.h @@ -29,13 +29,7 @@ * fraction of the sarcomere at sarcomere length @f$SL = SL_0 \, \lambda@f$ (with * @f$\lambda@f$ the fiber stretch), and @f$a_\text{XB}@f$ is the tension * upscaling factor. Because @f$\mu_P^1 + \mu_N^1@f$ and @f$\phi(SL)@f$ are - * dimensionless, @f$a_\text{XB}@f$ sets the units of the returned active tension; - * it is stored in the stress units of the simulation, so no separate output - * conversion is applied. - * - * The model is calibrated in a fixed unit system (calcium in [uM], time in [s], - * length in [um]); the svMultiPhysics inputs are converted to these units at the - * interface (see the conversion members below). + * dimensionless, @f$a_\text{XB}@f$ sets the units of the returned active tension. * * **References**: * 1. [Regazzoni, Dede', Quarteroni (2020)](https://doi.org/10.1371/journal.pcbi.1008294) @@ -82,37 +76,29 @@ class RegazzoniActiveStress : public ActiveStress { /** * @brief Model parameters class. * - * Declares the regulatory-unit (RU), crossbridge (XB), geometry, and tension - * parameters. The values registered below correspond to the published human - * body-temperature calibration of the reference implementation, expressed in - * the units documented for each member. The registered value of @c a_XB, - * 22.894, expresses the reference calibration in MPa; the value supplied in - * solver.xml must instead use the stress unit of the simulation's mechanical - * configuration. - * - * All parameters are required. A complete @c Regazzoni parameter block - * containing every parameter must be provided in solver.xml. Any value may be changed to - * use a different calibration, but omitting a parameter causes a parse error; - * the registered reference value is not used as an automatic default. + * 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 ActiveStressModelParameters { public: Parameters() : ActiveStressModelParameters(label) { constexpr bool required = true; - add_parameter("Kbasic", 13.0, required); - add_parameter("Koff", 100.0, required); + // Reference values: Regazzoni 2020 human body-temperature calibration, + // expressed consistently with the unit system used by this parameter set. + add_parameter("Kbasic", 0.013, required); + add_parameter("Koff", 0.1, required); add_parameter("Q", 2.0, required); add_parameter("mu", 10.0, required); add_parameter("gamma", 12.0, required); - add_parameter("Kd0", 0.381, required); - add_parameter("alphaKd", -0.571, required); + add_parameter("Kd0", 3.81e-4, required); + add_parameter("alphaKd", -5.71e-4, required); add_parameter("SL0", 2.2, required); - add_parameter("r0", 134.31, required); + add_parameter("r0", 0.13431, required); add_parameter("alpha", 25.184, required); - add_parameter("mu0_fP", 32.653, required); - add_parameter("mu1_fP", 0.778, required); + add_parameter("mu0_fP", 0.032653, required); + add_parameter("mu1_fP", 7.78e-4, required); add_parameter("LA", 1.25, required); add_parameter("LM", 1.65, required); @@ -226,7 +212,7 @@ class RegazzoniActiveStress : public ActiveStress { * effective boundary-neighbour transitions from the mean-field closure, then * updates @p state_RU in place. * - * @param[in] dt Substep size [s]. + * @param[in] dt Substep size [time]. * @param[in] rates_T Central-tropomyosin transition rates, * indexed @c rates_T[TL][TC][TR][CC]. * @param[in] rates_C Troponin transition rates, indexed @c rates_C[CC][TC]. @@ -245,8 +231,8 @@ class RegazzoniActiveStress : public ActiveStress { * transition rates from the updated RU probabilities, forms the 4x4 linear * system for the implicit update and solves it in place for @p state_XB. * - * @param[in] dt Outer time step [s]. - * @param[in] velocity Shortening velocity @f$-\dot{SL}/SL_0@f$ [s^-1]. + * @param[in] dt Outer time step [time]. + * @param[in] velocity Shortening velocity @f$-\dot{SL}/SL_0@f$ [1/time]. * @param[in] rates_T Central-tropomyosin transition rates, * indexed @c rates_T[TL][TC][TR][CC]. * @param[in] state_RU The updated 16 RU-state probabilities, @@ -266,61 +252,45 @@ class RegazzoniActiveStress : public ActiveStress { * thin and thick filaments overlap exactly once, a piecewise-linear function * of the sarcomere length built from the filament geometry (LA, LM, LB). * - * @param[in] sarcomere_length Sarcomere length @f$SL@f$ [um]. + * @param[in] sarcomere_length Sarcomere length @f$SL@f$ [length]. */ double fraction_single_overlap(double sarcomere_length) const; /// @} - /// @name Interface unit conversions (svMultiPhysics EM units to reference units) - /// @{ - - /// Calcium conversion, millimolar [mM] to micromolar [uM]. - static constexpr double calcium_mM_to_microM = 1.0e3; - - /// Time conversion, milliseconds [ms] to seconds [s]. - static constexpr double time_ms_to_s = 1.0e-3; - - /// @} - - /// RU forward-Euler substep size [s]. - static constexpr double ru_substep = 2.5e-5; + /// RU forward-Euler substep size [time]. + static constexpr double ru_substep = 2.5e-2; - /// Fixed reference sarcomere length [um] in the length-dependent dissociation + /// Fixed reference sarcomere length [length] in the length-dependent dissociation /// constant (distinct from the parameter SL0). static constexpr double kd_reference_sarcomere_length = 2.15; - /// @name RU model parameters (reference units) + /// @name RU model parameters /// @{ - double Kbasic; ///< Basic tropomyosin transition rate [s^-1]. - double Koff; ///< Troponin unbinding rate [s^-1]. + double Kbasic; ///< Basic tropomyosin transition rate [1/time]. + double Koff; ///< Troponin unbinding rate [1/time]. double Q; ///< Tropomyosin transition-rate asymmetry factor [-]. double mu; ///< Calcium-binding cooperativity factor [-]. double gamma; ///< Nearest-neighbour cooperativity factor [-]. - double Kd0; ///< Calcium dissociation constant at reference length [uM]. - double alphaKd; ///< Length dependence of the dissociation constant [uM/um]. - double SL0; ///< Reference sarcomere length [um]; maps stretch to length. + double Kd0; ///< Calcium dissociation constant at reference length [calcium]. + double alphaKd; ///< Length dependence of the dissociation constant [calcium/length]. + double SL0; ///< Reference sarcomere length [length]; maps stretch to length. - double r0; ///< Combined attachment-detachment rate at zero velocity [s^-1]. + double r0; ///< Combined attachment-detachment rate at zero velocity [1/time]. double alpha; ///< Coefficient of |v| in r(v) = r0 + alpha * |v| [-]. - double mu0_fP; ///< Permissive influx into the zeroth-moment crossbridge state [s^-1]. - double mu1_fP; ///< Permissive influx into the first-moment crossbridge state [s^-1]. + double mu0_fP; ///< Permissive influx into the zeroth-moment crossbridge state [1/time]. + double mu1_fP; ///< Permissive influx into the first-moment crossbridge state [1/time]. - double LA; ///< Thin-filament (actin) length [um]. - double LM; ///< Thick-filament (myosin) length [um]. - double LB; ///< Length of the myosin bare zone [um]. + double LA; ///< Thin-filament (actin) length [length]. + double LM; ///< Thick-filament (myosin) length [length]. + double LB; ///< Length of the myosin bare zone [length]. - /// Tension upscaling factor. + /// Tension upscaling factor [stress]. /// /// Because the crossbridge moments and the overlap fraction are dimensionless, - /// a_XB is the only quantity carrying stress units, so the returned active - /// tension has the same stress unit as a_XB and no stress-unit conversion is - /// performed. a_XB must therefore be expressed in the same stress unit as the - /// mechanical configuration. The reference calibration value 22.894 is - /// expressed in MPa, consistent with the coupled electromechanics slab case; - /// the equivalent values are 22.894e3 in kPa and 22.894e6 in Pa. Provide the - /// value matching the case's stress unit. + /// a_XB sets the stress unit of the returned active tension. It must be + /// expressed in the same stress unit as the mechanical configuration. double a_XB; /// @} diff --git a/tests/cases/electromechanics/slab_Regazzoni/result_001.vtu b/tests/cases/electromechanics/slab_Regazzoni/result_001.vtu index fc338238d..2817656f8 100644 --- a/tests/cases/electromechanics/slab_Regazzoni/result_001.vtu +++ b/tests/cases/electromechanics/slab_Regazzoni/result_001.vtu @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:31c6a24f61d33d7646970715f4e5c0dc359cd9b2fad077bdcfbd7625f428bbec -size 1443224 +oid sha256:6664a7689a0dc5c30f9b5238d8e3e0b9495ce0963cca4d525fd4200446647d52 +size 1444180 diff --git a/tests/cases/electromechanics/slab_Regazzoni/solver.xml b/tests/cases/electromechanics/slab_Regazzoni/solver.xml index ddc60007a..6d5d830e2 100644 --- a/tests/cases/electromechanics/slab_Regazzoni/solver.xml +++ b/tests/cases/electromechanics/slab_Regazzoni/solver.xml @@ -142,22 +142,22 @@ - 13.0 - 100.0 - 2.0 - 10.0 - 12.0 - 0.381 - -0.571 - 2.2 - 134.31 - 25.184 - 32.653 - 0.778 - 1.25 - 1.65 - 0.18 - 22.894 + 0.013 + 0.1 + 2.0 + 10.0 + 12.0 + 3.81e-4 + -5.71e-4 + 2.2 + 0.13431 + 25.184 + 0.032653 + 7.78e-4 + 1.25 + 1.65 + 0.18 + 22.894 From 4a9a0dec0ce93463ff689a976adf4ee63a5b2224 Mon Sep 17 00:00:00 2001 From: Sami Bismar Date: Wed, 5 Aug 2026 13:09:12 -0500 Subject: [PATCH 14/17] Revert Regazzoni active tension normalization --- .../Source/solver/active_stress_regazzoni.cpp | 13 ++------- Code/Source/solver/active_stress_regazzoni.h | 26 ++++++------------ .../electromechanics/slab_Regazzoni/README.md | 27 +++++++++---------- .../slab_Regazzoni/result_001.vtu | 4 +-- 4 files changed, 25 insertions(+), 45 deletions(-) diff --git a/Code/Source/solver/active_stress_regazzoni.cpp b/Code/Source/solver/active_stress_regazzoni.cpp index 3d8f2a78f..24c6fde0a 100644 --- a/Code/Source/solver/active_stress_regazzoni.cpp +++ b/Code/Source/solver/active_stress_regazzoni.cpp @@ -3,9 +3,7 @@ #include "active_stress_regazzoni.h" -#include "FE/Common/FEException.h" #include "eigen3/Eigen/Dense" -#include "utils.h" #include #include @@ -130,20 +128,13 @@ void RegazzoniActiveStress::advance_time_step_local( double RegazzoniActiveStress::compute_active_tension_local( const Vector &state, const double fiber_stretch) const { - if (utils::is_zero(fiber_stretch)) - svmp::raise( - "RegazzoniActiveStress: fiber_stretch is zero or near zero; this indicates a " - "degenerate (collapsed) element and is not a valid deformation state."); - const double sarcomere_length = SL0 * fiber_stretch; - // Compute paper tension T̃_act = a_XB * (μ_P^1 + μ_N^1) * φ(SL) from the + // Active tension T_act = a_XB * (μ_P^1 + μ_N^1) * φ(SL) from the // permissive and non-permissive XB first moments (state entries 17 and 19), // scaled by the single-overlap fraction and the upscaling factor a_XB. - // Return T̃_act / λ_f so that the assembled first-PK active stress has norm - // T̃_act, consistent with the RDQ20-MF paper formulation. return a_XB * (state[xb_index(1)] + state[xb_index(3)]) * - fraction_single_overlap(sarcomere_length) / fiber_stretch; + fraction_single_overlap(sarcomere_length); } void RegazzoniActiveStress::ru_transition_rates_tropomyosin( diff --git a/Code/Source/solver/active_stress_regazzoni.h b/Code/Source/solver/active_stress_regazzoni.h index be82e8c56..6d6fdb245 100644 --- a/Code/Source/solver/active_stress_regazzoni.h +++ b/Code/Source/solver/active_stress_regazzoni.h @@ -159,29 +159,19 @@ class RegazzoniActiveStress : public ActiveStress { Vector &state) const override; /** - * @brief Compute the PK2 active-tension coefficient for a single node. + * @brief Compute the scalar active tension for a single node. * - * Computes the paper tension - * @f$\widetilde{T}_\text{act} = a_\text{XB} (\mu_P^1 + \mu_N^1) \phi(SL)@f$ + * Returns the RDQ20-MF scalar active tension + * @f$T_\text{act} = a_\text{XB} (\mu_P^1 + \mu_N^1) \phi(SL)@f$ * from the XB first moments (state entries 17 and 19) and the single-overlap - * fraction at @f$SL = SL_0 \lambda_f@f$, then returns - * @f$\widetilde{T}_\text{act} / \lambda_f@f$. + * fraction at @f$SL = SL_0 \lambda_f@f$. Assembly of this scalar into the + * continuum active stress tensor follows the existing svMultiPhysics mechanics + * convention; the formulation of that assembly is outside the scope of this class. * - * The division by @f$\lambda_f@f$ is necessary because svMultiPhysics - * assembles the active stress as - * @f$\mathbf{S}_\text{act} = T_\text{code}\,\mathbf{f}_0\otimes\mathbf{f}_0@f$ - * (second-PK), so the resulting first-PK active stress satisfies - * @f$\|\mathbf{P}_\text{act}\| = \widetilde{T}_\text{act}@f$, matching the - * paper's definition of active tension. + * @p fiber_stretch is used to compute the sarcomere length + * @f$SL = SL_0 \lambda_f@f$ and therefore the overlap fraction @f$\phi(SL)@f$. * * The returned value has the stress units of @f$a_\text{XB}@f$. - * - * @note The slab regression test uses a (0.7, 0.2, 0.1) fiber/sheet/normal - * directional distribution, which is an svMultiPhysics extension of the - * paper's fiber-only active stress formulation. - * - * @throws svmp::FE::InvalidArgumentException if @p fiber_stretch is zero - * or near zero (degenerate/collapsed element). */ virtual double compute_active_tension_local(const Vector &state, diff --git a/tests/cases/electromechanics/slab_Regazzoni/README.md b/tests/cases/electromechanics/slab_Regazzoni/README.md index 57d8e92c4..bc1875964 100644 --- a/tests/cases/electromechanics/slab_Regazzoni/README.md +++ b/tests/cases/electromechanics/slab_Regazzoni/README.md @@ -56,25 +56,24 @@ are not prescribed by the RDQ20-MF model itself. ## Regression Reference -svMultiPhysics stores `T_code = T̃_act / λ_f` in each per-direction field -(`Active_tension_fibers`, `Active_tension_sheets`, `Active_tension_normal`), where -`T̃_act` is the paper tension from Regazzoni et al. (2020) and `λ_f` is the fiber -stretch; the paper tension is recovered as `T̃_act = λ_f · T_code`. Each -per-direction field stores `η · T_code`, where `η` is the corresponding directional -weight; their sum recovers `T_code`. +svMultiPhysics stores `T_act = a_XB * (μ_P^1 + μ_N^1) * φ(SL)` — the scalar +RDQ20-MF active tension — and distributes it into per-direction fields +(`Active_tension_fibers`, `Active_tension_sheets`, `Active_tension_normal`) using +the directional weights `η`. Each per-direction field stores `η · T_act`; their +sum recovers `T_act` because the directional weights sum to one. Assembly of +this scalar into the continuum active stress tensor follows the existing +svMultiPhysics mechanics convention. The formulation of that assembly will be +addressed separately. The active tension fields in `result_001.vtu` were validated node-by-node against the C++ reference implementation at commit [`26f05df`](https://github.com/FrancescoRegazzoni/cardiac-activation/commit/26f05df28891df7b3c69f16bb136cdced6b63c4d). Both implementations use the same implicit-Euler XB scheme, so agreement is -to machine precision (~1e-16 relative error). The comparison evaluates -`λ_f · T_code` (i.e. `T̃_act`) from the svMultiPhysics output against the reference -C++ active tension, using the calcium and sarcomere-length inputs from this one-step -test. Because this test covers only one time step from the undeformed configuration, -`λ_f ≈ 1.0` throughout, and the regression reference does not strongly distinguish -the normalized (`T̃_act / λ_f`) and unnormalized (`T̃_act`) conventions. -The remaining fields in the VTU serve as integrated svMultiPhysics regression -references and were not independently validated by the RDQ20-MF reference code. +to machine precision (~1e-16 relative error). The comparison evaluates `T_act` +from the svMultiPhysics output directly against the reference C++ active tension, +using the calcium and sarcomere-length inputs from this one-step test. The remaining +fields in the VTU serve as integrated svMultiPhysics regression references and were +not independently validated by the RDQ20-MF reference code. ## References diff --git a/tests/cases/electromechanics/slab_Regazzoni/result_001.vtu b/tests/cases/electromechanics/slab_Regazzoni/result_001.vtu index 2817656f8..c423af8b6 100644 --- a/tests/cases/electromechanics/slab_Regazzoni/result_001.vtu +++ b/tests/cases/electromechanics/slab_Regazzoni/result_001.vtu @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:6664a7689a0dc5c30f9b5238d8e3e0b9495ce0963cca4d525fd4200446647d52 -size 1444180 +oid sha256:d09218c7505190adf4963e2c080d4aeab5712d7f22376bfd53abc13388d7e7a4 +size 1441592 From 26667273e8b9765c23b95b2f7f776effb76a368e Mon Sep 17 00:00:00 2001 From: Sami Bismar Date: Wed, 5 Aug 2026 15:20:31 -0500 Subject: [PATCH 15/17] Make Regazzoni model parameters configurable --- .../Source/solver/active_stress_regazzoni.cpp | 4 ++ Code/Source/solver/active_stress_regazzoni.h | 37 ++++++++----------- .../slab_Regazzoni/solver.xml | 2 + 3 files changed, 21 insertions(+), 22 deletions(-) diff --git a/Code/Source/solver/active_stress_regazzoni.cpp b/Code/Source/solver/active_stress_regazzoni.cpp index 24c6fde0a..aeeb125e1 100644 --- a/Code/Source/solver/active_stress_regazzoni.cpp +++ b/Code/Source/solver/active_stress_regazzoni.cpp @@ -24,6 +24,8 @@ void RegazzoniActiveStress::read_model_specific_parameters( "and can produce a zero dissociation constant at physiological " "sarcomere lengths)."); SL0 = params.get_scalar("SL0"); + ru_substep = params.get_scalar("ru_substep"); + kd_reference_sarcomere_length = params.get_scalar("kd_reference_sarcomere_length"); r0 = params.get_scalar("r0"); alpha = params.get_scalar("alpha"); @@ -46,6 +48,8 @@ void RegazzoniActiveStress::distribute_model_specific_parameters( cm.bcast(cm_mod, &Kd0); cm.bcast(cm_mod, &alphaKd); cm.bcast(cm_mod, &SL0); + cm.bcast(cm_mod, &ru_substep); + cm.bcast(cm_mod, &kd_reference_sarcomere_length); cm.bcast(cm_mod, &r0); cm.bcast(cm_mod, &alpha); diff --git a/Code/Source/solver/active_stress_regazzoni.h b/Code/Source/solver/active_stress_regazzoni.h index 6d6fdb245..ef8260050 100644 --- a/Code/Source/solver/active_stress_regazzoni.h +++ b/Code/Source/solver/active_stress_regazzoni.h @@ -20,14 +20,14 @@ * moments with one implicit-Euler step per time step; the active tension is then * reconstructed from the XB first moments. * - * The active tension is + * The returned scalar active tension is * @f[ - * T_\text{act} = a_\text{XB} \, (\mu_P^1 + \mu_N^1) \, \phi(SL)\;, + * \Tact = a_\text{XB} \, (\mu_P^1 + \mu_N^1) \, \phi(SL)\;, * @f] * where @f$\mu_P^1@f$ and @f$\mu_N^1@f$ are the permissive and non-permissive * first XB moments (state entries 17 and 19), @f$\phi(SL)@f$ is the single-overlap - * fraction of the sarcomere at sarcomere length @f$SL = SL_0 \, \lambda@f$ (with - * @f$\lambda@f$ the fiber stretch), and @f$a_\text{XB}@f$ is the tension + * fraction of the sarcomere at sarcomere length @f$SL = SL_0 \, \fiberstretch@f$ + * (with @f$\fiberstretch@f$ the fiber stretch), and @f$a_\text{XB}@f$ is the tension * upscaling factor. Because @f$\mu_P^1 + \mu_N^1@f$ and @f$\phi(SL)@f$ are * dimensionless, @f$a_\text{XB}@f$ sets the units of the returned active tension. * @@ -94,6 +94,8 @@ class RegazzoniActiveStress : public ActiveStress { add_parameter("Kd0", 3.81e-4, required); add_parameter("alphaKd", -5.71e-4, required); add_parameter("SL0", 2.2, required); + add_parameter("ru_substep", 2.5e-2, required); + add_parameter("kd_reference_sarcomere_length", 2.15, required); add_parameter("r0", 0.13431, required); add_parameter("alpha", 25.184, required); @@ -161,17 +163,10 @@ class RegazzoniActiveStress : public ActiveStress { /** * @brief Compute the scalar active tension for a single node. * - * Returns the RDQ20-MF scalar active tension - * @f$T_\text{act} = a_\text{XB} (\mu_P^1 + \mu_N^1) \phi(SL)@f$ - * from the XB first moments (state entries 17 and 19) and the single-overlap - * fraction at @f$SL = SL_0 \lambda_f@f$. Assembly of this scalar into the - * continuum active stress tensor follows the existing svMultiPhysics mechanics - * convention; the formulation of that assembly is outside the scope of this class. - * - * @p fiber_stretch is used to compute the sarcomere length - * @f$SL = SL_0 \lambda_f@f$ and therefore the overlap fraction @f$\phi(SL)@f$. - * - * The returned value has the stress units of @f$a_\text{XB}@f$. + * Evaluates @f$\Tact@f$ as defined in the class description, using + * @p fiber_stretch to compute the sarcomere length + * @f$SL = SL_0 \, \fiberstretch@f$. The returned value has the stress + * units of @f$a_\text{XB}@f$. */ virtual double compute_active_tension_local(const Vector &state, @@ -248,13 +243,6 @@ class RegazzoniActiveStress : public ActiveStress { /// @} - /// RU forward-Euler substep size [time]. - static constexpr double ru_substep = 2.5e-2; - - /// Fixed reference sarcomere length [length] in the length-dependent dissociation - /// constant (distinct from the parameter SL0). - static constexpr double kd_reference_sarcomere_length = 2.15; - /// @name RU model parameters /// @{ @@ -266,6 +254,11 @@ class RegazzoniActiveStress : public ActiveStress { double Kd0; ///< Calcium dissociation constant at reference length [calcium]. double alphaKd; ///< Length dependence of the dissociation constant [calcium/length]. double SL0; ///< Reference sarcomere length [length]; maps stretch to length. + double ru_substep; ///< RU forward-Euler substep size [time]. + + /// Reference sarcomere length [length] used in the length-dependent + /// dissociation constant (distinct from the parameter SL0). + double kd_reference_sarcomere_length; double r0; ///< Combined attachment-detachment rate at zero velocity [1/time]. double alpha; ///< Coefficient of |v| in r(v) = r0 + alpha * |v| [-]. diff --git a/tests/cases/electromechanics/slab_Regazzoni/solver.xml b/tests/cases/electromechanics/slab_Regazzoni/solver.xml index 6d5d830e2..a55af9763 100644 --- a/tests/cases/electromechanics/slab_Regazzoni/solver.xml +++ b/tests/cases/electromechanics/slab_Regazzoni/solver.xml @@ -150,6 +150,8 @@ 3.81e-4 -5.71e-4 2.2 + 2.5e-2 + 2.15 0.13431 25.184 0.032653 From c289cf2432183cf6a393f67b1b0fbc6de783829d Mon Sep 17 00:00:00 2001 From: Sami Bismar Date: Wed, 5 Aug 2026 15:42:34 -0500 Subject: [PATCH 16/17] Replace Regazzoni C arrays with std::array --- .../Source/solver/active_stress_regazzoni.cpp | 45 +++++++++-------- Code/Source/solver/active_stress_regazzoni.h | 48 ++++++++++++------- 2 files changed, 55 insertions(+), 38 deletions(-) diff --git a/Code/Source/solver/active_stress_regazzoni.cpp b/Code/Source/solver/active_stress_regazzoni.cpp index aeeb125e1..fc376ff5d 100644 --- a/Code/Source/solver/active_stress_regazzoni.cpp +++ b/Code/Source/solver/active_stress_regazzoni.cpp @@ -76,8 +76,7 @@ void RegazzoniActiveStress::advance_time_step_local( const double sarcomere_length = SL0 * fiber_stretch; // Calcium/stretch-independent central-tropomyosin transition rates. - double rates_T[2][2][2][2]; - ru_transition_rates_tropomyosin(rates_T); + const RUArray rates_T = ru_transition_rates_tropomyosin(); // Troponin transition rates rates_C[CC][TC]: the calcium-binding row (CC = 0) // depends on calcium and sarcomere length; the unbinding row (CC = 1) does @@ -86,7 +85,7 @@ void RegazzoniActiveStress::advance_time_step_local( Koff / (Kd0 - alphaKd * (kd_reference_sarcomere_length - sarcomere_length)) * calcium; - double rates_C[2][2]; + BinaryPairArray rates_C; rates_C[0][0] = calcium_on_rate; rates_C[0][1] = calcium_on_rate; rates_C[1][0] = Koff; @@ -94,7 +93,7 @@ void RegazzoniActiveStress::advance_time_step_local( // Deserialize the 16 RU probabilities (entries 0-15). The crossbridge moments // (entries 16-19) are left untouched by this increment. - double state_RU[2][2][2][2]; + RUArray state_RU; for (int TL = 0; TL < 2; ++TL) for (int TC = 0; TC < 2; ++TC) for (int TR = 0; TR < 2; ++TR) @@ -113,10 +112,10 @@ void RegazzoniActiveStress::advance_time_step_local( // Advance the crossbridge moments (entries 16-19) from the updated RU state. // The velocity v = -dSL/dt / SL0 reduces to -d(lambda)/dt because SL = SL0 * lambda. const double velocity = -fiber_stretch_rate; - double state_XB[4]; + XBArray state_XB; for (int i = 0; i < 4; ++i) state_XB[i] = state[xb_index(i)]; - xb_implicit_update(dt, velocity, rates_T, state_RU, state_XB); + state_XB = xb_implicit_update(dt, velocity, rates_T, state_RU, state_XB); // Serialize the updated RU probabilities back into the state vector. for (int TL = 0; TL < 2; ++TL) @@ -141,8 +140,9 @@ double RegazzoniActiveStress::compute_active_tension_local( fraction_single_overlap(sarcomere_length); } -void RegazzoniActiveStress::ru_transition_rates_tropomyosin( - double (&rates_T)[2][2][2][2]) const { +RegazzoniActiveStress::RUArray +RegazzoniActiveStress::ru_transition_rates_tropomyosin() const { + RUArray rates_T; for (int TL = 0; TL < 2; ++TL) for (int TR = 0; TR < 2; ++TR) { const int permissive_neighbors = TL + TR; @@ -159,14 +159,15 @@ void RegazzoniActiveStress::ru_transition_rates_tropomyosin( rates_T[TL][0][TR][0] = opening_rate / mu; rates_T[TL][0][TR][1] = opening_rate; } + return rates_T; } void RegazzoniActiveStress::ru_forward_euler_substep( - double dt, const double (&rates_T)[2][2][2][2], - const double (&rates_C)[2][2], double (&state_RU)[2][2][2][2]) const { + double dt, const RUArray &rates_T, + const BinaryPairArray &rates_C, RUArray &state_RU) const { // Probability fluxes from central-unit transitions. - double flux_TC[2][2][2][2]; // central tropomyosin - double flux_CC[2][2][2][2]; // central troponin + RUArray flux_TC; // central tropomyosin + RUArray flux_CC; // central troponin for (int TL = 0; TL < 2; ++TL) for (int TC = 0; TC < 2; ++TC) for (int TR = 0; TR < 2; ++TR) @@ -180,8 +181,8 @@ void RegazzoniActiveStress::ru_forward_euler_substep( // Effective transition rates of the boundary neighbours, obtained from the // mean-field closure by conditioning the central-unit flux on the neighbour // pair state. - double rate_left[2][2]; - double rate_right[2][2]; + BinaryPairArray rate_left; + BinaryPairArray rate_right; for (int TL = 0; TL < 2; ++TL) for (int TC = 0; TC < 2; ++TC) { double flux_sum = 0.0; @@ -210,8 +211,8 @@ void RegazzoniActiveStress::ru_forward_euler_substep( // TR's only neighbour is TC on its left → rate_left[TC][TR]. // (rate_left == rate_right numerically due to mean-field LR symmetry, so the // result is unchanged, but the names now match the physical convention.) - double flux_TL[2][2][2][2]; // left tropomyosin - double flux_TR[2][2][2][2]; // right tropomyosin + RUArray flux_TL; // left tropomyosin + RUArray flux_TR; // right tropomyosin for (int TL = 0; TL < 2; ++TL) for (int TC = 0; TC < 2; ++TC) for (int TR = 0; TR < 2; ++TR) @@ -234,11 +235,11 @@ void RegazzoniActiveStress::ru_forward_euler_substep( flux_CC[TL][TC][TR][CC] + flux_CC[TL][TC][TR][1 - CC]); } -void RegazzoniActiveStress::xb_implicit_update( +RegazzoniActiveStress::XBArray RegazzoniActiveStress::xb_implicit_update( double dt, double velocity, - const double (&rates_T)[2][2][2][2], - const double (&state_RU)[2][2][2][2], - double (&state_XB)[4]) const { + const RUArray &rates_T, + const RUArray &state_RU, + const XBArray &state_XB) const { // Permissivity and the permissive/non-permissive probability fluxes from the // updated RU state. double permissivity = 0.0; @@ -290,8 +291,10 @@ void RegazzoniActiveStress::xb_implicit_update( const Eigen::Matrix solution = system.colPivHouseholderQr().solve(rhs); + XBArray result; for (int i = 0; i < 4; ++i) - state_XB[i] = solution(i); + result[i] = solution(i); + return result; } double RegazzoniActiveStress::fraction_single_overlap(double sarcomere_length) const { diff --git a/Code/Source/solver/active_stress_regazzoni.h b/Code/Source/solver/active_stress_regazzoni.h index ef8260050..426be475d 100644 --- a/Code/Source/solver/active_stress_regazzoni.h +++ b/Code/Source/solver/active_stress_regazzoni.h @@ -6,6 +6,8 @@ #include "active_stress.h" +#include + /** * @brief Mean-field active stress model (implements the RDQ20-MF formulation). * @@ -173,6 +175,16 @@ class RegazzoniActiveStress : public ActiveStress { const double fiber_stretch) const override; private: + /// Array indexed over the four binary RU configuration variables (TL, TC, TR, CC). + using RUArray = + std::array, 2>, 2>, 2>; + + /// Array indexed over a pair of binary state variables. + using BinaryPairArray = std::array, 2>; + + /// Array of the four crossbridge moment state variables. + using XBArray = std::array; + /// @name Regulatory-unit (RU) dynamics helpers /// @{ @@ -180,15 +192,15 @@ class RegazzoniActiveStress : public ActiveStress { * @brief Compute the central-tropomyosin transition rate for each local RU * configuration. * - * Fills @p rates_T, indexed as @c rates_T[TL][TC][TR][CC], where TL, TC and TR - * are the binary left, central and right tropomyosin states and CC is the - * central troponin calcium-binding state. Each entry is the rate at which the - * central tropomyosin changes state for that configuration. Because the rate - * depends on the neighbour states TL and TR, nearest-neighbour cooperativity - * is retained through the tracked TL-TC-TR configuration. These rates depend - * only on the model parameters, not on calcium or stretch. + * Returns an @c RUArray where entry @c [TL][TC][TR][CC] is the rate at which + * the central tropomyosin changes state for that configuration. Because the + * rate depends on the neighbour states TL and TR, nearest-neighbour + * cooperativity is retained through the tracked TL-TC-TR configuration. + * These rates depend only on the model parameters, not on calcium or stretch. + * + * @return Central-tropomyosin transition rates, indexed @c [TL][TC][TR][CC]. */ - void ru_transition_rates_tropomyosin(double (&rates_T)[2][2][2][2]) const; + RUArray ru_transition_rates_tropomyosin() const; /** * @brief Advance the 16 RU-state probabilities by one forward-Euler substep. @@ -205,16 +217,16 @@ class RegazzoniActiveStress : public ActiveStress { * indexed @c state_RU[TL][TC][TR][CC]. */ void ru_forward_euler_substep(double dt, - const double (&rates_T)[2][2][2][2], - const double (&rates_C)[2][2], - double (&state_RU)[2][2][2][2]) const; + const RUArray &rates_T, + const BinaryPairArray &rates_C, + RUArray &state_RU) const; /** * @brief Advance the four crossbridge moments by one implicit-Euler step. * * Computes the permissivity and the effective permissive/non-permissive * transition rates from the updated RU probabilities, forms the 4x4 linear - * system for the implicit update and solves it in place for @p state_XB. + * system for the implicit update, and returns the updated moments. * * @param[in] dt Outer time step [time]. * @param[in] velocity Shortening velocity @f$-\dot{SL}/SL_0@f$ [1/time]. @@ -222,13 +234,15 @@ class RegazzoniActiveStress : public ActiveStress { * indexed @c rates_T[TL][TC][TR][CC]. * @param[in] state_RU The updated 16 RU-state probabilities, * indexed @c state_RU[TL][TC][TR][CC]. - * @param[in,out] state_XB The four crossbridge moments, ordered + * @param[in] state_XB The four crossbridge moments (input), ordered + * @f$[\mu_P^0, \mu_P^1, \mu_N^0, \mu_N^1]@f$. + * @return Updated crossbridge moments, ordered * @f$[\mu_P^0, \mu_P^1, \mu_N^0, \mu_N^1]@f$. */ - void xb_implicit_update(double dt, double velocity, - const double (&rates_T)[2][2][2][2], - const double (&state_RU)[2][2][2][2], - double (&state_XB)[4]) const; + XBArray xb_implicit_update(double dt, double velocity, + const RUArray &rates_T, + const RUArray &state_RU, + const XBArray &state_XB) const; /** * @brief Single-overlap fraction of the sarcomere at a given length. From cb56643e0c9d837f3bb6d44e99cf84d8005028a7 Mon Sep 17 00:00:00 2001 From: Sami Bismar Date: Wed, 5 Aug 2026 15:55:08 -0500 Subject: [PATCH 17/17] Consolidate electromechanics slab test configurations --- tests/cases/electromechanics/slab/README.md | 94 +++++++++++++++--- ...lt_001.vtu => result_NashPanfilov_001.vtu} | 0 .../result_Regazzoni_001.vtu} | 0 .../{solver.xml => solver_NashPanfilov.xml} | 0 .../solver.xml => slab/solver_Regazzoni.xml} | 0 .../electromechanics/slab_Regazzoni/README.md | 98 ------------------- .../slab_Regazzoni/mesh/X0.vtp | 3 - .../slab_Regazzoni/mesh/X1.vtp | 3 - .../slab_Regazzoni/mesh/volume.vtu | 3 - tests/test_electromechanics.py | 12 +-- 10 files changed, 85 insertions(+), 128 deletions(-) rename tests/cases/electromechanics/slab/{result_001.vtu => result_NashPanfilov_001.vtu} (100%) rename tests/cases/electromechanics/{slab_Regazzoni/result_001.vtu => slab/result_Regazzoni_001.vtu} (100%) rename tests/cases/electromechanics/slab/{solver.xml => solver_NashPanfilov.xml} (100%) rename tests/cases/electromechanics/{slab_Regazzoni/solver.xml => slab/solver_Regazzoni.xml} (100%) delete mode 100644 tests/cases/electromechanics/slab_Regazzoni/README.md delete mode 100644 tests/cases/electromechanics/slab_Regazzoni/mesh/X0.vtp delete mode 100644 tests/cases/electromechanics/slab_Regazzoni/mesh/X1.vtp delete mode 100644 tests/cases/electromechanics/slab_Regazzoni/mesh/volume.vtu diff --git a/tests/cases/electromechanics/slab/README.md b/tests/cases/electromechanics/slab/README.md index c37f66d62..383876bc6 100755 --- a/tests/cases/electromechanics/slab/README.md +++ b/tests/cases/electromechanics/slab/README.md @@ -1,13 +1,23 @@ # **Problem Description** -Simulate cardiac electromechanics on a slab of myocardial tissue. This test -couples 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. +Simulate cardiac electromechanics on a slab of myocardial tissue. This +directory contains two solver configurations that share the same geometry and +electrophysiology setup but differ in the active-stress model: -## Electrophysiology +| Configuration file | Active-stress model | +|------------------------------|---------------------| +| `solver_NashPanfilov.xml` | Nash-Panfilov | +| `solver_Regazzoni.xml` | RDQ20-MF (Regazzoni)| + +Both 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. + +## Shared Geometry and Electrophysiology + +The mesh is a rectangular slab (`mesh/`) with two boundary faces `X0` and `X1`. The propagation of the transmembrane potential is modeled with the ten-Tusscher-Panfilov (`TTP`) cell activation model [2, 3], using epicardial @@ -25,13 +35,16 @@ into two `Domain`s: an unstimulated region (`domain 1`) and a stimulated region ``` -## Mechanics - The tissue is modeled as a nearly incompressible Holzapfel-Ogden material with -modified anisotropy (`HolzapfelOgden-ModifiedAnisotropy`) [4]. Active contraction -is driven by the calcium concentration computed by the electrophysiology model, -through the Nash-Panfilov active-stress model [5] with a directional distribution -along the fiber, sheet, and sheet-normal directions. +modified anisotropy (`HolzapfelOgden-ModifiedAnisotropy`) [4]. The slab is fixed +with a zero-displacement Dirichlet boundary condition on the `X1` face, and +contracts as the depolarization wave propagates through the tissue. + +## Nash-Panfilov variant (`solver_NashPanfilov.xml`) + +Active contraction is driven by the calcium concentration computed by the +electrophysiology model, through the Nash-Panfilov active-stress model [5] with a +directional distribution along the fiber, sheet, and sheet-normal directions. ``` @@ -45,8 +58,55 @@ along the fiber, sheet, and sheet-normal directions. ``` -The slab is fixed with a zero-displacement Dirichlet boundary condition on the -`X1` face, and contracts as the depolarization wave propagates through the tissue. +**Regression reference:** `result_NashPanfilov_001.vtu` + +## Regazzoni variant (`solver_Regazzoni.xml`) + +Active contraction is driven by the calcium concentration computed by the +electrophysiology model, through the RDQ20-MF mean-field active-stress model [6], +configured with the published human body-temperature calibration expressed in the +solver's unit system (time in ms, calcium in mM, length in µm). The scalar active +tension is distributed along the fiber, sheet, and sheet-normal directions using +the same directional weights as the Nash-Panfilov variant. + +``` + + Regazzoni + + 0.7 + 0.2 + 0.1 + + ... + +``` + +The (0.7, 0.2, 0.1) fiber/sheet/sheet-normal directional weights are an +svMultiPhysics extension of the paper's fiber-only active stress formulation; they +are not prescribed by the RDQ20-MF model itself. + +**Regression reference:** `result_Regazzoni_001.vtu` + +### Validation + +svMultiPhysics stores `T_act = a_XB * (μ_P^1 + μ_N^1) * φ(SL)` — the scalar +RDQ20-MF active tension — and distributes it into per-direction fields +(`Active_tension_fibers`, `Active_tension_sheets`, `Active_tension_normal`) using +the directional weights `η`. Each per-direction field stores `η · T_act`; their +sum recovers `T_act` because the directional weights sum to one. Assembly of +this scalar into the continuum active stress tensor follows the existing +svMultiPhysics mechanics convention. The formulation of that assembly will be +addressed separately. + +The active tension fields in `result_Regazzoni_001.vtu` were validated +node-by-node against the C++ reference implementation at commit +[`26f05df`](https://github.com/FrancescoRegazzoni/cardiac-activation/commit/26f05df28891df7b3c69f16bb136cdced6b63c4d). +Both implementations use the same implicit-Euler XB scheme, so agreement is +to machine precision (~1e-16 relative error). The comparison evaluates `T_act` +from the svMultiPhysics output directly against the reference C++ active tension, +using the calcium and sarcomere-length inputs from this one-step test. The remaining +fields in the VTU serve as integrated svMultiPhysics regression references and were +not independently validated by the RDQ20-MF reference code. ## References @@ -68,4 +128,8 @@ of the Royal Society A, 367(1902):3445–3475, 2009. [5] M. P. Nash and A. V. Panfilov. Electromechanical model of excitable tissue to study reentrant cardiac arrhythmias. Progress in Biophysics and Molecular Biology, -85(2-3):501–522, 2004. \ No newline at end of file +85(2-3):501–522, 2004. + +[6] F. Regazzoni, L. Dede', and A. Quarteroni. Biophysically detailed mathematical +models of multiscale cardiac active mechanics. PLOS Computational Biology, +16(10):e1008294, 2020. diff --git a/tests/cases/electromechanics/slab/result_001.vtu b/tests/cases/electromechanics/slab/result_NashPanfilov_001.vtu similarity index 100% rename from tests/cases/electromechanics/slab/result_001.vtu rename to tests/cases/electromechanics/slab/result_NashPanfilov_001.vtu diff --git a/tests/cases/electromechanics/slab_Regazzoni/result_001.vtu b/tests/cases/electromechanics/slab/result_Regazzoni_001.vtu similarity index 100% rename from tests/cases/electromechanics/slab_Regazzoni/result_001.vtu rename to tests/cases/electromechanics/slab/result_Regazzoni_001.vtu diff --git a/tests/cases/electromechanics/slab/solver.xml b/tests/cases/electromechanics/slab/solver_NashPanfilov.xml similarity index 100% rename from tests/cases/electromechanics/slab/solver.xml rename to tests/cases/electromechanics/slab/solver_NashPanfilov.xml diff --git a/tests/cases/electromechanics/slab_Regazzoni/solver.xml b/tests/cases/electromechanics/slab/solver_Regazzoni.xml similarity index 100% rename from tests/cases/electromechanics/slab_Regazzoni/solver.xml rename to tests/cases/electromechanics/slab/solver_Regazzoni.xml diff --git a/tests/cases/electromechanics/slab_Regazzoni/README.md b/tests/cases/electromechanics/slab_Regazzoni/README.md deleted file mode 100644 index bc1875964..000000000 --- a/tests/cases/electromechanics/slab_Regazzoni/README.md +++ /dev/null @@ -1,98 +0,0 @@ - -# **Problem Description** - -Simulate cardiac electromechanics on a slab of myocardial tissue. This test -couples 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. - -## Electrophysiology - -The propagation of the transmembrane potential is modeled with the -ten-Tusscher-Panfilov (`TTP`) cell activation model [2, 3], using epicardial -parameters (included from `../../cep/ttp_parameters/ttp_epicardium_parameters.xml`) -and anisotropic conductivity aligned with the fiber direction. The domain is split -into two `Domain`s: an unstimulated region (`domain 1`) and a stimulated region -(`domain 2`) where an external `Istim` stimulus initiates depolarization. - -``` - - -35.714 - 0.0 - 2.0 - 10000.0 - -``` - -## Mechanics - -The tissue is modeled as a nearly incompressible Holzapfel-Ogden material with -modified anisotropy (`HolzapfelOgden-ModifiedAnisotropy`) [4]. Active contraction -is driven by the calcium concentration computed by the electrophysiology model, -through the RDQ20-MF mean-field active-stress model [5], configured with the -published human body-temperature calibration. The scalar active tension is -distributed along the fiber, sheet, and sheet-normal directions using the same -directional weights as the NashPanfilov slab case. - -``` - - Regazzoni - - 0.7 - 0.2 - 0.1 - - ... - -``` - -The slab is fixed with a zero-displacement Dirichlet boundary condition on the -`X1` face, and contracts as the depolarization wave propagates through the tissue. - -The (0.7, 0.2, 0.1) fiber/sheet/sheet-normal directional weights are an -svMultiPhysics extension of the paper's fiber-only active stress formulation; they -are not prescribed by the RDQ20-MF model itself. - -## Regression Reference - -svMultiPhysics stores `T_act = a_XB * (μ_P^1 + μ_N^1) * φ(SL)` — the scalar -RDQ20-MF active tension — and distributes it into per-direction fields -(`Active_tension_fibers`, `Active_tension_sheets`, `Active_tension_normal`) using -the directional weights `η`. Each per-direction field stores `η · T_act`; their -sum recovers `T_act` because the directional weights sum to one. Assembly of -this scalar into the continuum active stress tensor follows the existing -svMultiPhysics mechanics convention. The formulation of that assembly will be -addressed separately. - -The active tension fields in `result_001.vtu` were validated node-by-node against -the C++ reference implementation at commit -[`26f05df`](https://github.com/FrancescoRegazzoni/cardiac-activation/commit/26f05df28891df7b3c69f16bb136cdced6b63c4d). -Both implementations use the same implicit-Euler XB scheme, so agreement is -to machine precision (~1e-16 relative error). The comparison evaluates `T_act` -from the svMultiPhysics output directly against the reference C++ active tension, -using the calcium and sarcomere-length inputs from this one-step test. The remaining -fields in the VTU serve as integrated svMultiPhysics regression references and were -not independently validated by the RDQ20-MF reference code. - -## References - -[1] S. A. Niederer, E. Kerfoot, A. P. Benson, et al. Verification of cardiac tissue -electrophysiology simulators using an N-version benchmark. Philosophical Transactions -of the Royal Society A, 369(1954):4331–4351, 2011. - -[2] K. H. W. J. ten Tusscher, D. Noble, P. J. Noble, and A. V. Panfilov. A model for -human ventricular tissue. American Journal of Physiology-Heart and Circulatory -Physiology, 286(4):H1573–H1589, apr 2004. - -[3] K. H. W. J. ten Tusscher and A. V. Panfilov. Alternans and spiral breakup in a -human ventricular tissue model. American Journal of Physiology-Heart and Circulatory -Physiology, 291(3):H1088–H1100, sep 2006. - -[4] G. A. Holzapfel and R. W. Ogden. Constitutive modelling of passive myocardium: a -structurally based framework for material characterization. Philosophical Transactions -of the Royal Society A, 367(1902):3445–3475, 2009. - -[5] F. Regazzoni, L. Dede', and A. Quarteroni. Biophysically detailed mathematical -models of multiscale cardiac active mechanics. PLOS Computational Biology, -16(10):e1008294, 2020. diff --git a/tests/cases/electromechanics/slab_Regazzoni/mesh/X0.vtp b/tests/cases/electromechanics/slab_Regazzoni/mesh/X0.vtp deleted file mode 100644 index eaebcc0eb..000000000 --- a/tests/cases/electromechanics/slab_Regazzoni/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_Regazzoni/mesh/X1.vtp b/tests/cases/electromechanics/slab_Regazzoni/mesh/X1.vtp deleted file mode 100644 index 07ad82435..000000000 --- a/tests/cases/electromechanics/slab_Regazzoni/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_Regazzoni/mesh/volume.vtu b/tests/cases/electromechanics/slab_Regazzoni/mesh/volume.vtu deleted file mode 100644 index dd12aca89..000000000 --- a/tests/cases/electromechanics/slab_Regazzoni/mesh/volume.vtu +++ /dev/null @@ -1,3 +0,0 @@ -version https://git-lfs.github.com/spec/v1 -oid sha256:2d9c504f6d9dc221739d107e19414233cdda640d73647eba03902a1605274149 -size 421967 diff --git a/tests/test_electromechanics.py b/tests/test_electromechanics.py index 92a5d8bb0..c728e3185 100644 --- a/tests/test_electromechanics.py +++ b/tests/test_electromechanics.py @@ -1,6 +1,4 @@ from .conftest import run_with_reference -import os -import subprocess # Common folder for all tests in this file base_folder = "electromechanics" @@ -24,10 +22,12 @@ def test_slab(n_proc): - test_folder = "slab" - run_with_reference(base_folder, test_folder, fields, n_proc, t_max=1) + run_with_reference(base_folder, "slab", fields, n_proc, t_max=1, + name_inp="solver_NashPanfilov.xml", + name_ref="result_NashPanfilov_001.vtu") def test_slab_regazzoni(n_proc): - test_folder = "slab_Regazzoni" - run_with_reference(base_folder, test_folder, fields, n_proc, t_max=1) + run_with_reference(base_folder, "slab", fields, n_proc, t_max=1, + name_inp="solver_Regazzoni.xml", + name_ref="result_Regazzoni_001.vtu")