Add Land Niederer Active Stress Model for Electromechanics

Code/Source/solver/ActiveStressLandNiederer.cpp

@@ -0,0 +1,156 @@
+// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the
+// University of California, and others. SPDX-License-Identifier: BSD-3-Clause
+
+#include "ActiveStressLandNiederer.h"
+
+void ActiveStressLandNiederer::read_model_specific_parameters(
+    const ActiveStressModelParameters &params) {
+  ActiveStressODE::read_model_specific_parameters(params);
+
+  calcium_scaling_factor = params.get_scalar("calcium_scaling_factor");
+  CaRef = params.get_scalar("CaRef");
+  eta_Tm = params.get_scalar("eta_Tm");
+  k_uw = params.get_scalar("k_uw");
+  k_ws = params.get_scalar("k_ws");
+  Tref = params.get_scalar("Tref");
+  k_TRPN = params.get_scalar("k_TRPN");
+  eta_TRPN = params.get_scalar("eta_TRPN");
+  k_u = params.get_scalar("k_u");
+  TRPN50 = params.get_scalar("TRPN50");
+  rw = params.get_scalar("rw");
+  rs = params.get_scalar("rs");
+  gamma_s = params.get_scalar("gamma_s");
+  gamma_w = params.get_scalar("gamma_w");
+  phi = params.get_scalar("phi");
+  Aeff = params.get_scalar("Aeff");
+  beta0 = params.get_scalar("beta0");
+  beta1 = params.get_scalar("beta1");
+
+  disable_force_strain_rate_feedback_ =
+      params.get_bool("Disable_force_strain_rate_feedback");
+
+  // With the feedback disabled the stretch rate computation can be skipped.
+  needs_fiber_stretch_rate_ = !disable_force_strain_rate_feedback_;
+}
+
+void ActiveStressLandNiederer::distribute_model_specific_parameters(const CmMod &cm_mod,
+                                                        const cmType &cm) {
+  ActiveStressODE::distribute_model_specific_parameters(cm_mod, cm);
+
+  cm.bcast(cm_mod, &calcium_scaling_factor);
+  cm.bcast(cm_mod, &CaRef);
+  cm.bcast(cm_mod, &eta_Tm);
+  cm.bcast(cm_mod, &k_uw);
+  cm.bcast(cm_mod, &k_ws);
+  cm.bcast(cm_mod, &Tref);
+  cm.bcast(cm_mod, &k_TRPN);
+  cm.bcast(cm_mod, &eta_TRPN);
+  cm.bcast(cm_mod, &k_u);
+  cm.bcast(cm_mod, &TRPN50);
+  cm.bcast(cm_mod, &rw);
+  cm.bcast(cm_mod, &rs);
+  cm.bcast(cm_mod, &gamma_s);
+  cm.bcast(cm_mod, &gamma_w);
+  cm.bcast(cm_mod, &phi);
+  cm.bcast(cm_mod, &Aeff);
+  cm.bcast(cm_mod, &beta0);
+  cm.bcast(cm_mod, &beta1);
+  cm.bcast(cm_mod, &disable_force_strain_rate_feedback_);
+  cm.bcast(cm_mod, &needs_fiber_stretch_rate_);
+}
+
+void ActiveStressLandNiederer::init_local(Vector<double> &state) const { 
+    state[0] = 1.0; 
+    state[1] = 0.0; 
+    state[2] = 0.0; 
+    state[3] = 0.0; 
+    state[4] = 0.0; 
+    state[5] = 0.0; 
+    const double lambda = 1.0; 
+    state[6] = std::max(0.0, 1.0 + beta0 * (lambda + std::min(0.87, lambda) - 1.87)); 
+}
+
+void ActiveStressLandNiederer::advance_time_step_local(const double t, const double dt,
+                                           const double calcium,
+                                           const double fiber_stretch,
+                                           const double fiber_stretch_rate,
+                                           Vector<double> &state) const {
+
+  Vector<double> ode_state(6); 
+  for (int i = 0; i < 6; i++) {
+    ode_state[i] = state[i]; 
+  }
+
+  ActiveStressODE::advance_time_step_local(t, dt, calcium, fiber_stretch,
+                                           fiber_stretch_rate, ode_state);
+
+  for (unsigned int i = 0; i < 6; ++i) {
+    state[i] = ode_state[i];
+  }
+
+  const double lambda = std::min(1.2, fiber_stretch); 
+  state[6] = std::max(0.0, 1.0 + beta0 * (lambda + std::min(0.87, lambda) - 1.87));
+}
+
+Vector<double> ActiveStressLandNiederer::getf(const double t, const Vector<double> &state,
+                                  const double calcium,
+                                  const double fiber_stretch,
+                                  const double fiber_stretch_rate) const {
+  Vector<double> f(6);
+
+  // State Variables
+  const double XB = state[0]; 
+  const double XW = std::max(0.0, state[1]); 
+  const double CaTRPN = std::max(0.0, state[2]); 
+  const double XS = std::max(0.0, state[3]); 
+  const double ZETAW = state[4]; 
+  const double ZETAS = state[5]; 
+
+  // Bound State Equation
+  const double lambda = std::min(1.2, fiber_stretch); 
+  const double XU = (1 - XB) - XS - XW; 
+  const double k_b = k_u * std::pow(TRPN50, eta_Tm)/ (1 - rs - (1 - rs)*rw);
+  f[0] = k_b * std::min(100.0, std::pow(CaTRPN, -(eta_Tm/2.0))) * XU  - k_u * std::pow(CaTRPN, (eta_Tm/2.0)) *  XB;
+
+  // Pre-power Stroke Kinetics 
+  const double k_wu = (k_uw * (1.0/rw - 1.0) - k_ws);
+  const double gamma_wu = gamma_w * std::abs(ZETAW); 
+  f[1] = k_uw*XU - k_wu*XW - k_ws*XW - gamma_wu*XW;
+
+  // Calcium-Troponin Binding Kinetics  
+  const double CaT50 = CaRef + beta1*std::min(0.2, lambda - 1.0); 
+  f[2] = k_TRPN * (std::pow(((calcium*calcium_scaling_factor)/CaT50), eta_TRPN) * (1.0 - CaTRPN) - CaTRPN); 
+
+  // Post-power Stroke Kinetics                                 
+  const double k_su = k_ws * rw * (1.0/rs - 1.0); 
+  const double term1 = (ZETAS > 0.0) ? ZETAS : 0.0;
+  const double term2 = (ZETAS < -1.0) ? (-ZETAS - 1.0) : 0.0;
+  const double gamma_su = gamma_s * std::max(term1, term2);
+  f[3] = k_ws*XW - k_su*XS - gamma_su*XS; 
+
+  // Distortion-Decay Kinetics
+  const double stretch_rate =
+      disable_force_strain_rate_feedback_ ? 0.0 : fiber_stretch_rate;
+  const double Aw = Aeff * rs/((1.0 - rs) * rw + rs);
+  const double cw = phi * k_uw * (1.0 - rs)*(1.0 - rw) / ((1.0 - rs)*rw);
+  f[4] = Aw * stretch_rate - cw*ZETAW;
+
+  const double As = Aw;
+  const double cs = phi * k_ws * (1.0 - rs)*rw/rs;
+  f[5] =  As * stretch_rate - cs*ZETAS;
+
+  return f;
+}
+
+double ActiveStressLandNiederer::compute_active_tension_local(const Vector<double> &state,
+                                                   const double fiber_stretch) const {
+  const double XW = std::max(0.0, state[1]);
+  const double XS = std::max(0.0, state[3]); 
+  const double ZETAW = state[4];  
+  const double ZETAS = state[5];
+  const double LFac = state[6];  
+
+  return LFac * (Tref/rs) * ((ZETAS+1.0) * XS + (ZETAW) * XW);
+}
+
+REGISTER_ACTIVE_STRESS_MODEL("LandNiederer", ActiveStressLandNiederer);

Code/Source/solver/ActiveStressLandNiederer.h

@@ -0,0 +1,223 @@
+// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the
+// University of California, and others. SPDX-License-Identifier: BSD-3-Clause
+
+#ifndef ACTIVE_STRESS_LAND_NIEDERER_H
+#define ACTIVE_STRESS_LAND_NIEDERER_H
+
+#include "ActiveStressODE.h"
+
+/**
+ * @brief Land-Niederer active stress model.
+ *
+ * This class implements the phenomenological Land-Niederer active stress model 
+ * [1], which represents the cross-bridge cycling and calcium-troponin binding 
+ * kinetics in a single set of ODEs. The cross-bridges are modeled as a three-state 
+ * system, with transitions between unbound (U), pre-power stroke (W), and post-power
+ * stroke (S) states. The blocked state (B) represents the tropomyosin blocking 
+ * the actin binding sites, and the unblocked state (U) repesents the availability
+ * of these sites for cross-bridge binding. The unblocking of the tropomyosin is 
+ * regulated by the calcium-troponin binding state (CaTRPN) and the fiber stretch. 
+ * The distortion of the cross-bridges are also represented for the pre-power stroke
+ * and post-power stroke states. The active tension is computed as a function of the 
+ * cross-bridge states and their distortions, scaled by a reference tension Tref. 
+ * This is given as
+ * @f[
+ *  \Tact = \frac{Tref}{r_s} \left[ (1 + \zeta_S) X_S + \zeta_W X_W \right]\;,
+ * @f]
+ * where @f$X_S@f$ and @f$X_W@f$ are the fractions of cross-bridges in the post-power 
+ * stroke and pre-power stroke states, respectively. @f$\zeta_S@f$ and @f$\zeta_W@f$ 
+ * are the distortions, and @f$r_s@f$ is the steady-state duty ratio of the 
+ * cross-bridges. 
+ * **References**:
+ * 1. [Land, Niederer (2017)](https://doi.org/10.1016/j.yjmcc.2017.03.008)
+ */
+class ActiveStressLandNiederer : public ActiveStressODE {
+public:
+  /// Model label, used for factory registration and XML selection. 
+  static inline const std::string label = "LandNiederer";
+
+  /*
+   * @brief Model parameters class.
+   * Declares the parameters required by the model. All parameters are 
+   * marked as required, and omitting a parameter will cause a parse error.
+  */ 
+  class Parameters : public ActiveStressODE::Parameters {
+  public:
+    Parameters() : ActiveStressODE::Parameters(label) {
+      constexpr bool required = true;
+
+      // Reference values calibrated in Land 2017 
+      add_parameter("calcium_scaling_factor", 1000.0, required);
+      add_parameter("CaRef", 0.805, required);
+      add_parameter("eta_Tm", 5.0, required);
+      add_parameter("k_uw", 0.182, required);
+      add_parameter("k_ws", 0.012, required);
+      add_parameter("Tref", 120.0, required);
+      add_parameter("k_TRPN", 0.1, required);
+      add_parameter("eta_TRPN", 2.0, required);
+      add_parameter("k_u", 1.0, required); 
+      add_parameter("TRPN50", 0.35, required);
+      add_parameter("rw", 0.5, required);
+      add_parameter("rs", 0.25, required);
+      add_parameter("gamma_s", 0.0085, required);
+      add_parameter("gamma_w", 0.615, required);
+      add_parameter("phi", 2.23, required);
+      add_parameter("Aeff", 25.0, required);
+      add_parameter("beta0", 2.3, required);
+      add_parameter("beta1", -2.4, required);
+
+      add_parameter("Disable_force_strain_rate_feedback", false, !required);
+    }
+  };
+
+  /**
+   * @brief Constructor.
+   */
+  ActiveStressLandNiederer() : ActiveStressODE(/* n_state_variables = */ 7,
+                                   /* needs_fiber_stretch = */ true,
+                                   /* needs_fiber_stretch_rate = */ true) {}
+
+  /**
+   * @brief Construct an instance of model parameters.
+   */
+  virtual std::unique_ptr<ActiveStressModelParameters>
+  get_parameters() const override {
+    return std::make_unique<Parameters>();
+  }
+
+protected:
+  /**
+   * @brief Read model parameters from a parameter object.
+   */
+  virtual void read_model_specific_parameters(
+      const ActiveStressModelParameters &params) override;
+
+  /**
+   * @brief Distribute model parameters to all parallel processes.
+   */
+  virtual void distribute_model_specific_parameters(const CmMod &cm_mod,
+                                                    const cmType &cm) override;
+
+  /**
+   * @brief Initialize the state vector for a single node.
+   *
+   * @param[out] state State vector for a single node, to be initialized by
+   *   this function.
+   */
+  virtual void init_local(Vector<double> &state) const override;
+
+  /**
+   * @brief Advance in time for a single node.
+   *
+   */
+  virtual void advance_time_step_local(const double t, const double dt,
+                                       const double calcium,
+                                       const double fiber_stretch,
+                                       const double fiber_stretch_rate,
+                                       Vector<double> &state) const override;
+
+  /**
+   * @brief Compute the rate of change in the state variables.
+   */
+  virtual Vector<double> getf(const double t, const Vector<double> &state,
+                              const double calcium, const double fiber_stretch,
+                              const double fiber_stretch_rate) const override;
+
+  /**
+   * @brief Compute the active tension for a single node.
+   */
+  virtual double
+  compute_active_tension_local(const Vector<double> &state,
+                               const double fiber_stretch) const override;
+
+  /// @name Model parameters.
+  /// @{
+  /// Scaling factor required if converting calcium concentration from 
+  /// ionic model to contraction model 
+  double calcium_scaling_factor; 
+
+  /// Reference intracellular calcium concentration giving half-maximal
+  /// troponin C saturation, @f$[Ca^{2+}]_{T50,ref}@f$, used in the CaTRPN 
+  /// binding ODE [uM]
+  double CaRef;
+
+  /// Cooperativity (Hill) exponent @f$n_{Tm}@f$ for the tropomyosin
+  /// blocked/unblocked (B/U) transition; sets the steepness of thin-filament
+  /// activation by CaTRPN [-]
+  double eta_Tm;
+
+  /// Rate constant for transition from unbound (U) to pre-power stroke (W) 
+  /// state [1/ms]
+  double k_uw;
+
+  /// Rate constant for transition from pre-power stroke (W) to post-power 
+  /// stroke (S) state [1/ms]
+  double k_ws;
+
+  /// Maximum observable tension @f$\Tref@f$ at resting length [kPa]
+  double Tref;
+
+  /// Rate constant @f$k_{TRPN}@f$ governing calcium binding/unbinding
+  /// kinetics of troponin C [1/ms]
+  double k_TRPN;
+
+  /// Cooperativity (Hill) exponent @f$n_{TRPN}@f$ of the calcium-troponin C
+  /// binding rate [-]
+  double eta_TRPN;
+
+  /// Tropomyosin unblocking rate constant @f$k_u@f$ (rate at which blocked
+  /// binding sites become unblocked); together with rw, rs, and TRPN50 it
+  /// fixes the blocking rate kb [1/ms]
+  double k_u;
+
+  /// The value of CaTRPN where bounded state B = 0.5 in steady-state [-]
+  double TRPN50;
+
+  /// Steady-state ratio @f$r_w@f$ between the pre-powerstroke (W) population
+  /// and the non-strongly-bound (U+W) population at equilibrium; sets the
+  /// reverse rate constant kwu [-]
+  double rw;
+
+  /// Steady-state duty ratio @f$r_s@f$: the fraction of cross-bridges in the
+  /// strongly-bound (post-powerstroke, S) state at equilibrium; sets the
+  /// reverse rate constant ksu [-]
+  double rs;
+
+  /// Strain-dependent detachment-rate coefficient @f$\gamma_s@f$ for the
+  /// strongly-bound (post-powerstroke) S state; scales how fast cross-bridges
+  /// are pulled off by distortion in that state [-]
+  double gamma_s;
+
+  /// Strain-dependent detachment-rate coefficient @f$\gamma_w@f$ for the
+  /// weakly-bound (pre-powerstroke) W state; scales how fast cross-bridges
+  /// are pulled off by distortion in that state [-]
+  double gamma_w;
+
+  /// Proportionality factor @f$\phi@f$ relating the cross-bridge distortion
+  /// decay rates (cw, cs) to the forward cycling rates k_uw, k_ws [-]
+  double phi;
+
+  /// Rescaled magnitude @f$A_{eff}@f$ of the immediate cross-bridge
+  /// distortion response to fiber stretch rate; sets the distortion
+  /// sensitivities As, Aw [-]
+  double Aeff;
+
+  /// Coefficient @f$\beta_0@f$ controlling the length-dependence of maximal
+  /// tension through changes in filament overlap (Frank-Starling effect) [-]
+  double beta0;
+
+  /// Coefficient @f$\beta_1@f$ controlling the length-dependence of calcium
+  /// sensitivity (shifts CaRef/[Ca2+]T50 with sarcomere stretch) [-]
+  double beta1;
+
+  /// Controls force-strain-rate feedback in the cross-bridge distortion ODEs.
+  /// - @c false (default): the ZETAW/ZETAS distortion states respond to the
+  ///   fiber stretch rate as in Land, Niederer (2017).
+  /// - @c true: the fiber stretch rate contribution is zeroed out, disabling
+  ///   the feedback.
+  bool disable_force_strain_rate_feedback_ = false;
+
+  /// @}
+};
+
+#endif