diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt index c5ab81146..b37e8b985 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_regazzoni.cpp + SPLIT.c svZeroD_interface/LPNSolverInterface.h svZeroD_interface/LPNSolverInterface.cpp 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_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.cpp b/Code/Source/solver/active_stress_regazzoni.cpp new file mode 100644 index 000000000..fc376ff5d --- /dev/null +++ b/Code/Source/solver/active_stress_regazzoni.cpp @@ -0,0 +1,315 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause + +#include "active_stress_regazzoni.h" + +#include "eigen3/Eigen/Dense" + +#include +#include + +void RegazzoniActiveStress::read_model_specific_parameters( + const ActiveStressModelParameters ¶ms) { + 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"); + if (alphaKd > 0.0) + svmp::raise( + "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)."); + 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"); + 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 RegazzoniActiveStress::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); + 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); + cm.bcast(cm_mod, &ru_substep); + cm.bcast(cm_mod, &kd_reference_sarcomere_length); + + cm.bcast(cm_mod, &r0); + 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 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 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 { + const double sarcomere_length = SL0 * fiber_stretch; + + // Calcium/stretch-independent central-tropomyosin transition rates. + 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 + // not. + const double calcium_on_rate = + Koff / + (Kd0 - alphaKd * (kd_reference_sarcomere_length - sarcomere_length)) * + calcium; + BinaryPairArray rates_C; + 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. + RUArray state_RU; + 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 - 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 velocity v = -dSL/dt / SL0 reduces to -d(lambda)/dt because SL = SL0 * lambda. + const double velocity = -fiber_stretch_rate; + XBArray state_XB; + for (int i = 0; i < 4; ++i) + state_XB[i] = state[xb_index(i)]; + 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) + 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 RegazzoniActiveStress::compute_active_tension_local( + const Vector &state, const double fiber_stretch) const { + const double sarcomere_length = SL0 * fiber_stretch; + + // 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 a_XB * (state[xb_index(1)] + state[xb_index(3)]) * + fraction_single_overlap(sarcomere_length); +} + +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; + + // 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; + } + return rates_T; +} + +void RegazzoniActiveStress::ru_forward_euler_substep( + double dt, const RUArray &rates_T, + const BinaryPairArray &rates_C, RUArray &state_RU) const { + // Probability fluxes from central-unit transitions. + 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) + 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. + 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; + 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. + // 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.) + 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) + for (int CC = 0; CC < 2; ++CC) { + flux_TL[TL][TC][TR][CC] = + state_RU[TL][TC][TR][CC] * rate_right[TC][TL]; + flux_TR[TL][TC][TR][CC] = + state_RU[TL][TC][TR][CC] * rate_left[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]); +} + +RegazzoniActiveStress::XBArray RegazzoniActiveStress::xb_implicit_update( + double dt, double velocity, + 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; + 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); + XBArray result; + for (int i = 0; i < 4; ++i) + result[i] = solution(i); + return result; +} + +double RegazzoniActiveStress::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("Regazzoni", RegazzoniActiveStress); diff --git a/Code/Source/solver/active_stress_regazzoni.h b/Code/Source/solver/active_stress_regazzoni.h new file mode 100644 index 000000000..426be475d --- /dev/null +++ b/Code/Source/solver/active_stress_regazzoni.h @@ -0,0 +1,296 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause + +#ifndef ACTIVE_STRESS_REGAZZONI_H +#define ACTIVE_STRESS_REGAZZONI_H + +#include "active_stress.h" + +#include + +/** + * @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, 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 — 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 returned scalar active tension is + * @f[ + * \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 \, \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. + * + * **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 RegazzoniActiveStress : public ActiveStress { +public: + /// Model label, used for factory registration and XML selection. + static inline const std::string label = "Regazzoni"; + + /// @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 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 unsigned int xb_index(unsigned int i) { + return n_ru_states + i; + } + + /// @} + + /** + * @brief Model parameters class. + * + * Declares the parameters required by the model. All parameters are + * marked as required, and omitting a parameter will cause a parse error. + */ + class Parameters : public ActiveStressModelParameters { + public: + Parameters() : ActiveStressModelParameters(label) { + constexpr bool required = true; + + // 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", 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); + 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); + add_parameter("LB", 0.18, required); + add_parameter("a_XB", 22.894, required); + } + }; + + /** + * @brief Constructor. + */ + RegazzoniActiveStress() : 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. + * + * 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, + const double fiber_stretch, + const double fiber_stretch_rate, + Vector &state) const override; + + /** + * @brief Compute the scalar active tension for a single node. + * + * 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, + 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 + /// @{ + + /** + * @brief Compute the central-tropomyosin transition rate for each local RU + * configuration. + * + * 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]. + */ + RUArray ru_transition_rates_tropomyosin() 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 [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]. + * @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 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 returns the updated moments. + * + * @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, + * indexed @c state_RU[TL][TC][TR][CC]. + * @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$. + */ + 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. + * + * 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$ [length]. + */ + double fraction_single_overlap(double sarcomere_length) const; + + /// @} + + /// @name RU model parameters + /// @{ + + 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 [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| [-]. + 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 [length]. + double LM; ///< Thick-filament (myosin) length [length]. + double LB; ///< Length of the myosin bare zone [length]. + + /// Tension upscaling factor [stress]. + /// + /// Because the crossbridge moments and the overlap fraction are dimensionless, + /// 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; + + /// @} +}; + +#endif 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; 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/result_Regazzoni_001.vtu b/tests/cases/electromechanics/slab/result_Regazzoni_001.vtu new file mode 100644 index 000000000..c423af8b6 --- /dev/null +++ b/tests/cases/electromechanics/slab/result_Regazzoni_001.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:d09218c7505190adf4963e2c080d4aeab5712d7f22376bfd53abc13388d7e7a4 +size 1441592 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/solver_Regazzoni.xml b/tests/cases/electromechanics/slab/solver_Regazzoni.xml new file mode 100644 index 000000000..a55af9763 --- /dev/null +++ b/tests/cases/electromechanics/slab/solver_Regazzoni.xml @@ -0,0 +1,197 @@ + + + + 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 + + + + Regazzoni + + + 0.7 + 0.2 + 0.1 + + + + 0.013 + 0.1 + 2.0 + 10.0 + 12.0 + 3.81e-4 + -5.71e-4 + 2.2 + 2.5e-2 + 2.15 + 0.13431 + 25.184 + 0.032653 + 7.78e-4 + 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..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,5 +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): + run_with_reference(base_folder, "slab", fields, n_proc, t_max=1, + name_inp="solver_Regazzoni.xml", + name_ref="result_Regazzoni_001.vtu") 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