diff --git a/Code/Source/solver/ActiveStressLandNiederer.cpp b/Code/Source/solver/ActiveStressLandNiederer.cpp
new file mode 100644
index 000000000..b4295df4b
--- /dev/null
+++ b/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);
\ No newline at end of file
diff --git a/Code/Source/solver/ActiveStressLandNiederer.h b/Code/Source/solver/ActiveStressLandNiederer.h
new file mode 100644
index 000000000..82448315b
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/Code/Source/solver/ActiveStressODE.cpp b/Code/Source/solver/ActiveStressODE.cpp
index 9b9ca9e17..6fa2b92b1 100644
--- a/Code/Source/solver/ActiveStressODE.cpp
+++ b/Code/Source/solver/ActiveStressODE.cpp
@@ -9,6 +9,8 @@ void ActiveStressODE::read_model_specific_parameters(
 
   if (solver_str == "FE") {
     ode_solver = ODESolver::ForwardEuler;
+  } else if (solver_str == "RK4" || solver_str == "RK") {
+    ode_solver = ODESolver::RungeKutta4;
   } else {
     svmp::raise<svmp::ParseException>("Unknown ODE solver " + solver_str +
                                       " for active stress models.");
@@ -29,6 +31,19 @@ void ActiveStressODE::advance_time_step_local(const double t, const double dt,
     const Vector<double> f =
         getf(t - dt, state, calcium, fiber_stretch, fiber_stretch_rate);
     state.add(dt, f);
+  } else if (ode_solver == ODESolver::RungeKutta4) {
+    const Vector<double> k1 =
+        getf(t - dt, state, calcium, fiber_stretch, fiber_stretch_rate);
+    const Vector<double> k2 =
+        getf(t - 0.5 * dt, state + 0.5 * dt * k1, calcium, fiber_stretch,
+             fiber_stretch_rate);
+    const Vector<double> k3 =
+        getf(t - 0.5 * dt, state + 0.5 * dt * k2, calcium, fiber_stretch,
+             fiber_stretch_rate);
+    const Vector<double> k4 = getf(t, state + dt * k3, calcium, fiber_stretch,
+                                   fiber_stretch_rate);
+
+    state.add(dt / 6.0, k1 + 2.0 * (k2 + k3) + k4);
   } else {
     svmp::raise<svmp::InternalErrorException>(
         "Unknown ODE solver " + std::to_string(static_cast<int>(ode_solver)) +
diff --git a/Code/Source/solver/ActiveStressODE.h b/Code/Source/solver/ActiveStressODE.h
index 3efd47a31..6b573d71e 100644
--- a/Code/Source/solver/ActiveStressODE.h
+++ b/Code/Source/solver/ActiveStressODE.h
@@ -69,7 +69,34 @@ class ActiveStressODE : public ActiveStress {
      *                                    \fiberstretchrate_i^n)\;.
      * @f]
      */
-    ForwardEuler
+    ForwardEuler,
+
+    /**
+     * @brief 4th order explicit Runge-Kutta.
+     *
+     * The state vector is updated as follows:
+     * @f[ \begin{aligned}
+     *  \mathbf{k}_1 &= \mathbf{F}_\text{AS}(t^n, \astressstate_i^n,
+     *    \calcium_i^n, \fiberstretch_i^n, \fiberstretchrate_i^n)\;, \\
+     *  \mathbf{k}_2 &= \mathbf{F}_\text{AS}(t^n + \Delta t / 2,
+     *    \astressstate_i^n + \Delta t \mathbf{k}_1 / 2, \calcium_i^n,
+     *    \fiberstretch_i^n, \fiberstretchrate_i^n)\;, \\
+     *  \mathbf{k}_3 &= \mathbf{F}_\text{AS}(t^n + \Delta t / 2,
+     *    \astressstate_i^n + \Delta t \mathbf{k}_2 / 2, \calcium_i^n,
+     *    \fiberstretch_i^n, \fiberstretchrate_i^n)\;, \\
+     *  \mathbf{k}_4 &= \mathbf{F}_\text{AS}(t^{n+1},
+     *    \astressstate_i^n + \Delta t \mathbf{k}_3, \calcium_i^n,
+     *    \fiberstretch_i^n, \fiberstretchrate_i^n)\;, \\
+     *  \astressstate_i^{n+1} &= \astressstate_i^n + \frac{\Delta t}{6}
+     *    \left( \mathbf{k}_1 + 2 \mathbf{k}_2 + 2 \mathbf{k}_3
+     *           + \mathbf{k}_4 \right)\;.
+     * \end{aligned} @f]
+     *
+     * Note that the driving quantities @f$\calcium@f$, @f$\fiberstretch@f$ and
+     * @f$\fiberstretchrate@f$ are held fixed at their values at @f$t^{n+1}@f$
+     * across all four stages.
+     */
+    RungeKutta4
   };
 
   /**
@@ -125,11 +152,9 @@ class ActiveStressODE : public ActiveStress {
    * @param[in,out] state State vector for a single node, to be updated by
    *   this function.
    *
-   * @todo[michelebucelli] It might be necessary or useful to implement other
-   *   timestepping schemes, e.g. Runge-Kutta. In that case, we might want to
-   *   expand the interface to support implicit time stepping too, e.g. by
-   *   adding a method to evaluate the Jacobian matrix of the system, as in
-   *   @ref IonicModel.
+   * @todo[michelebucelli] It might be necessary or useful to implement
+   *   implicit time stepping too, e.g. by adding a method to evaluate the
+   *   Jacobian matrix of the system, as in @ref IonicModel.
    */
   virtual void advance_time_step_local(const double t, const double dt,
                                        const double calcium,
diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt
index 3f57b4daf..9e1968895 100644
--- a/Code/Source/solver/CMakeLists.txt
+++ b/Code/Source/solver/CMakeLists.txt
@@ -270,6 +270,7 @@ set(CSRCS
   ActiveStressODE.cpp
   ActiveStressNashPanfilov.cpp
   ActiveStressRegazzoni.cpp
+  ActiveStressLandNiederer.cpp
 
   SPLIT.c
 
diff --git a/tests/cases/electromechanics/slab/README.md b/tests/cases/electromechanics/slab/README.md
index 383876bc6..d5167435a 100755
--- a/tests/cases/electromechanics/slab/README.md
+++ b/tests/cases/electromechanics/slab/README.md
@@ -2,15 +2,16 @@
 # **Problem Description**
 
 Simulate cardiac electromechanics on a slab of myocardial tissue. This
-directory contains two solver configurations that share the same geometry and
+directory contains three solver configurations that share the same geometry and
 electrophysiology setup but differ in the active-stress model:
 
 | Configuration file           | Active-stress model |
 |------------------------------|---------------------|
 | `solver_NashPanfilov.xml`    | Nash-Panfilov       |
 | `solver_Regazzoni.xml`       | RDQ20-MF (Regazzoni)|
+| `solver_LandNiederer.xml`    | Land-Niederer       | 
 
-Both configurations couple cardiac electrophysiology (`CEP`) to solid mechanics
+The configurations couple cardiac electrophysiology (`CEP`) to solid mechanics
 (`struct`), reproducing the geometry and stimulation setting of the Niederer
 electrophysiology benchmark [1] with the addition of active contraction and
 finite-strain mechanics.
@@ -87,6 +88,28 @@ are not prescribed by the RDQ20-MF model itself.
 
 **Regression reference:** `result_Regazzoni_001.vtu`
 
+## Land-Niederer variant (`solver_LandNiederer.xml`) 
+
+Active contraction is driven by the calcium concentration computed by the
+electrophysiology model, through the Land-Niederer active-stress model [7]. The 
+model parameters are those calibrated to human ventricular cardiomyocyte data. The
+scalar active tension is distributed along the fiber, sheet, and sheet-normal 
+directions using the same directional weights as the Nash-Panfilov variant.
+
+```
+<Active_stress>
+  <Model>LandNiederer</Model>
+  <Directional_distribution>
+    <Fiber_direction> 0.7 </Fiber_direction>
+    <Sheet_direction> 0.2 </Sheet_direction>
+    <Sheet_normal_direction> 0.1 </Sheet_normal_direction>
+  </Directional_distribution>
+  ...
+</Active_stress>
+```
+
+**Regression reference:** `result_LandNiederer_001.vtu`
+
 ### Validation
 
 svMultiPhysics stores `T_act = a_XB * (μ_P^1 + μ_N^1) * φ(SL)` — the scalar
@@ -108,6 +131,11 @@ using the calcium and sarcomere-length inputs from this one-step test. The remai
 fields in the VTU serve as integrated svMultiPhysics regression references and were
 not independently validated by the RDQ20-MF reference code.
 
+The active tension fields in  `result_LandNiederer_001.vtu` were validated
+node-by-node against the MATLAB reference implementation provided 
+[here](https://www.cemrg.co.uk/models). A custom script with a python wrapper 
+for the MATLAB function is available in this [repository](https://github.com/kko27/generate_ref_solution.git).  
+
 ## References
 
 [1] S. A. Niederer, E. Kerfoot, A. P. Benson, et al. Verification of cardiac tissue
@@ -133,3 +161,7 @@ study reentrant cardiac arrhythmias. Progress in Biophysics and Molecular Biolog
 [6] F. Regazzoni, L. Dede', and A. Quarteroni. Biophysically detailed mathematical
 models of multiscale cardiac active mechanics. PLOS Computational Biology,
 16(10):e1008294, 2020.
+
+[7] S. Land, S. Park-Holohan, N. P. Smith, et al. A model of cardiac contraction
+based on novel measurements of tension development in human cardiomyocytes. 
+Journal of Molecular and Cellular Cardiology, 106:68-83, apr 2017. 
\ No newline at end of file
diff --git a/tests/cases/electromechanics/slab/result_LandNiederer_001.vtu b/tests/cases/electromechanics/slab/result_LandNiederer_001.vtu
new file mode 100644
index 000000000..1af9d63e5
--- /dev/null
+++ b/tests/cases/electromechanics/slab/result_LandNiederer_001.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:804a6ecf2b5bacb1a8e55392a6302138652e58f883a17ee10c99e4397cad334d
+size 1470461
diff --git a/tests/cases/electromechanics/slab/solver_LandNiederer.xml b/tests/cases/electromechanics/slab/solver_LandNiederer.xml
new file mode 100644
index 000000000..5e57517d3
--- /dev/null
+++ b/tests/cases/electromechanics/slab/solver_LandNiederer.xml
@@ -0,0 +1,204 @@
+<?xml version="1.0" encoding="UTF-8" ?>
+<svMultiPhysicsFile version="0.1">
+  <GeneralSimulationParameters>
+    <Continue_previous_simulation> false </Continue_previous_simulation>
+    <Number_of_spatial_dimensions> 3 </Number_of_spatial_dimensions>
+    <Number_of_time_steps> 1 </Number_of_time_steps>
+    <Time_step_size> 1.0 </Time_step_size>
+    <Spectral_radius_of_infinite_time_step> 0.50 </Spectral_radius_of_infinite_time_step>
+    <Searched_file_name_to_trigger_stop> STOP_SIM </Searched_file_name_to_trigger_stop>
+
+    <Save_results_to_VTK_format> true </Save_results_to_VTK_format>
+    <Name_prefix_of_saved_VTK_files> result </Name_prefix_of_saved_VTK_files>
+    <Increment_in_saving_VTK_files> 1 </Increment_in_saving_VTK_files>
+    <Start_saving_after_time_step> 0 </Start_saving_after_time_step>
+
+    <Increment_in_saving_restart_files> 1000 </Increment_in_saving_restart_files>
+    <Convert_BIN_to_VTK_format> 0 </Convert_BIN_to_VTK_format>
+
+    <Verbose> 1 </Verbose>
+    <Warning> 1 </Warning>
+    <Debug> 0 </Debug>
+  </GeneralSimulationParameters>
+
+  <Add_mesh name="block">
+    <Mesh_file_path> ./mesh/volume.vtu </Mesh_file_path>
+
+    <Add_face name="X0">
+      <Face_file_path> ./mesh/X0.vtp </Face_file_path>
+    </Add_face>
+
+    <Add_face name="X1">
+      <Face_file_path> ./mesh/X1.vtp </Face_file_path>
+    </Add_face>
+
+    <Domain_file_path> ./mesh/volume.vtu </Domain_file_path>
+
+    <Fiber_direction> (1, 0, 0)  </Fiber_direction>
+    <Fiber_direction> (0, 1, 0)  </Fiber_direction>
+    <Fiber_direction> (0, 1, 0)  </Fiber_direction>
+  </Add_mesh>
+
+  <Add_equation type="CEP">
+    <Coupled> true </Coupled>
+
+    <Min_iterations> 1 </Min_iterations>
+    <Max_iterations> 1 </Max_iterations>
+    <Tolerance> 1e-12 </Tolerance>
+
+    <Domain id="1">
+      <Electrophysiology_model> TTP </Electrophysiology_model>
+
+      <Isotropic_conductivity> 0.012571 </Isotropic_conductivity>
+      <Anisotropic_conductivity> 0.082715 </Anisotropic_conductivity>
+      <Anisotropic_conductivity> 0.0 </Anisotropic_conductivity>
+      <Anisotropic_conductivity> 0.0 </Anisotropic_conductivity>
+
+      <Include_xml> ../../cep/ttp_parameters/ttp_epicardium_parameters.xml </Include_xml>
+
+      <TTP>
+        <G_Na> 14.838 </G_Na>
+        <G_CaL> 3.98E-5 </G_CaL>
+        <G_Kr> 0.153 </G_Kr>
+      </TTP>
+
+      <ODE_solver> RK4 </ODE_solver>
+    </Domain>
+
+    <Domain id="2">
+      <Electrophysiology_model> TTP </Electrophysiology_model>
+
+      <Isotropic_conductivity> 0.012571 </Isotropic_conductivity>
+      <Anisotropic_conductivity> 0.082715 </Anisotropic_conductivity>
+      <Anisotropic_conductivity> 0.0 </Anisotropic_conductivity>
+      <Anisotropic_conductivity> 0.0 </Anisotropic_conductivity>
+
+      <Include_xml> ../../cep/ttp_parameters/ttp_epicardium_parameters.xml </Include_xml>
+
+      <TTP>
+        <G_Na> 14.838 </G_Na>
+        <G_CaL> 3.98E-5 </G_CaL>
+        <G_Kr> 0.153 </G_Kr>
+      </TTP>
+
+      <Stimulus type="Istim" >
+        <Amplitude> -35.714 </Amplitude>
+        <Start_time> 0.0 </Start_time>
+        <Duration> 2.0 </Duration>
+        <Cycle_length> 10000.0 </Cycle_length>
+      </Stimulus>
+
+      <ODE_solver> RK4 </ODE_solver>
+    </Domain>
+
+    <Output type="Spatial" >
+      <Membrane_potential> true </Membrane_potential>
+      <Ionic_state_variables> true </Ionic_state_variables>
+    </Output>
+
+    <LS type="CG" >
+      <Linear_algebra type="fsils" >
+        <Preconditioner> fsils </Preconditioner>
+      </Linear_algebra>
+      <Max_iterations> 100 </Max_iterations>
+      <Tolerance> 1e-12 </Tolerance>
+      <Krylov_space_dimension> 50 </Krylov_space_dimension>
+    </LS>
+  </Add_equation>
+
+  <Add_equation type="struct">
+    <Min_iterations> 1 </Min_iterations>
+    <Max_iterations> 6 </Max_iterations>
+    <Tolerance> 1e-12 </Tolerance>
+
+    <Density> 1e-3 </Density> <!-- g/mm^3 -->
+
+    <Constitutive_model type="HolzapfelOgden-ModifiedAnisotropy">
+      <a> 59.0e-6 </a>     <!-- MPa -->
+      <b> 8.023 </b>
+      <a4f> 18472.0e-6 </a4f> <!-- MPa -->
+      <b4f> 16.026 </b4f>
+      <a4s> 2481.0e-6 </a4s> <!-- MPa -->
+      <b4s> 11.12 </b4s>
+      <afs> 216.0e-6 </afs>  <!-- MPa -->
+      <bfs> 11.436 </bfs>
+      <k> 100.0 </k>
+    </Constitutive_model>
+
+    <Dilational_penalty_model> ST91 </Dilational_penalty_model>
+    <Penalty_parameter> 1.0 </Penalty_parameter> <!-- MPa -->
+
+    <Viscosity model="Potential" >
+        <Value> 1.0 </Value>
+    </Viscosity>
+
+    <Active_stress>
+      <Model>LandNiederer</Model>
+
+      <Directional_distribution>
+        <Fiber_direction> 0.7 </Fiber_direction>
+        <Sheet_direction> 0.2 </Sheet_direction>
+        <Sheet_normal_direction> 0.1 </Sheet_normal_direction>
+      </Directional_distribution>
+
+      <LandNiederer>
+        <ODE_solver> RK4 </ODE_solver>
+
+        <!-- 1000 as conversion factor from mM (TTP-EP model) to μM -->
+        <calcium_scaling_factor> 1000.0 </calcium_scaling_factor>   
+        <CaRef> 0.805 </CaRef>    <!-- μM -->
+        <eta_Tm> 5.0 </eta_Tm>    <!-- [-] -->
+        <k_uw> 0.182 </k_uw>      <!-- 1/ms -->
+        <k_ws> 0.012 </k_ws>      <!-- 1/ms -->
+        <Tref> 0.12 </Tref>       <!-- MPa -->
+        <k_TRPN> 0.1 </k_TRPN>    <!-- 1/ms -->
+        <eta_TRPN> 2.0 </eta_TRPN> <!-- [-] -->
+        <k_u> 1.0 </k_u>          <!-- 1/ms -->
+        <TRPN50> 0.35 </TRPN50>   <!-- [-] -->
+        <rw> 0.5 </rw>            <!-- [-] -->
+        <rs> 0.25 </rs>           <!-- [-] -->
+        <gamma_s> 0.0085 </gamma_s> <!-- [-] -->
+        <gamma_w> 0.615 </gamma_w> <!-- [-] -->
+        <phi> 2.23 </phi>          <!-- [-] -->
+        <Aeff> 25.0 </Aeff>        <!-- [-] -->
+        <beta0> 2.3 </beta0>       <!-- [-] -->
+        <beta1> -2.4 </beta1>      <!-- μM -->
+
+        <Disable_force_strain_rate_feedback> true </Disable_force_strain_rate_feedback>
+      </LandNiederer>
+    </Active_stress>
+
+    <Output type="Spatial" >
+      <Displacement> true </Displacement>
+      <Velocity> true </Velocity>
+      <Jacobian> true </Jacobian>
+      <Stress> true </Stress>
+      <Strain> true </Strain>
+      <Cauchy_stress> true </Cauchy_stress>
+      <Def_grad> true </Def_grad>
+      <VonMises_stress> true </VonMises_stress>
+
+      <Active_tension_fibers> true </Active_tension_fibers>
+      <Active_tension_sheets> true </Active_tension_sheets>
+      <Active_tension_normal> true </Active_tension_normal>
+
+      <Fiber_stretch> true </Fiber_stretch>
+    </Output>
+
+    <LS type="GMRES" >
+      <Linear_algebra type="fsils" >
+        <Preconditioner> fsils </Preconditioner>
+      </Linear_algebra>
+      <Tolerance> 1e-12 </Tolerance>
+      <Absolute_tolerance> 1e-14 </Absolute_tolerance>
+      <Max_iterations> 1000 </Max_iterations>
+    </LS>
+
+    <Add_BC name="X1" >
+      <Type> Dir </Type>
+      <Value> 0.0 </Value>
+    </Add_BC>
+  </Add_equation>
+</svMultiPhysicsFile>
+
+
diff --git a/tests/test_electromechanics.py b/tests/test_electromechanics.py
index e4f63bbd0..6f79c1749 100644
--- a/tests/test_electromechanics.py
+++ b/tests/test_electromechanics.py
@@ -23,7 +23,7 @@
 ]
 
 
-@pytest.mark.parametrize("model", ["NashPanfilov", "Regazzoni"])
+@pytest.mark.parametrize("model", ["NashPanfilov", "Regazzoni", "LandNiederer"])
 def test_slab(model, n_proc):
     run_with_reference(base_folder, "slab", fields, n_proc, t_max=1,
                        name_inp=f"solver_{model}.xml",
diff --git a/tests/unitTests/active_stress_tests/test_active_stress_land_niederer.cpp b/tests/unitTests/active_stress_tests/test_active_stress_land_niederer.cpp
new file mode 100644
index 000000000..a4cc64e3e
--- /dev/null
+++ b/tests/unitTests/active_stress_tests/test_active_stress_land_niederer.cpp
@@ -0,0 +1,38 @@
+// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the
+// University of California, and others. SPDX-License-Identifier: BSD-3-Clause
+
+/// @file
+
+#include "ActiveStressLandNiederer.h"
+#include "active_stress_test_helpers.h"
+#include "gtest/gtest.h"
+
+/**
+ * @test Run a standalone ActiveStressLandNiederer twitch for 600 updates with
+ * @f$\Delta t=1\,\mathrm{ms}@f$, the prescribed calcium and fiber-stretch
+ * trajectories, and the multiscale Land 2017 calibration configured below
+ * (svMP's model defaults).
+ *
+ * The trusted reference is generated by driving the original MATLAB
+ * @c Land.m implementation (adapted from https://www.cemrg.co.uk/models) 
+ * through the MATLAB engine with a manual Forward-Euler integrator that 
+ * mirrors @c ActiveStressODE's update rule; it is not generated by svMultiPhysics. 
+ * See this [repo](https://github.com/kko27/generate_ref_solution.git) to access the 
+ * reference solution generator and model. See @ref ActiveStressLandNiederer 
+ * for the model equations and state ordering: states 0-5 are the Land 2017 ODE states 
+ * (XB, XW, CaTRPN, XS, ZETAW, ZETAS) and state 6 is the deterministic Frank-Starling 
+ * length-dependence factor LFac.
+ */
+TEST(ActiveStressTrajectory, LandNiederer) {
+  ActiveStressLandNiederer::Parameters params;
+
+  ActiveStressTrajectoryConfiguration configuration;
+  configuration.final_time = 600.0;
+  configuration.time_step = 1.0;
+  configuration.reference_csv_filename =
+      "active_stress_land_niederer_twitch.csv";
+
+  ActiveStressTrajectoryTest<ActiveStressLandNiederer> trajectory(
+      params, configuration);
+  trajectory.run();
+}
diff --git a/tests/unitTests/reference_data/active_stress_land_niederer_twitch.csv b/tests/unitTests/reference_data/active_stress_land_niederer_twitch.csv
new file mode 100644
index 000000000..ac7a17f61
--- /dev/null
+++ b/tests/unitTests/reference_data/active_stress_land_niederer_twitch.csv
@@ -0,0 +1,12 @@
+step,Ta,s0,s1,s2,s3,s4,s5,s6
+0,0.0000000000000000e+00,1.0000000000000000e+00,0.0000000000000000e+00,1.5431503414220131e-03,0.0000000000000000e+00,0.0000000000000000e+00,0.0000000000000000e+00,1.0000000000000000e+00
+10,3.8835726655655308e-05,9.9999002345980692e-01,2.9941229181942714e-06,1.0517193396842484e-02,8.0907763865948556e-08,0.0000000000000000e+00,0.0000000000000000e+00,1.0000000000000000e+00
+30,2.0300434401505334e+00,5.9239113723449799e-01,1.2882747028414332e-01,4.5983992372518667e-01,4.2360968011456853e-03,-5.1401253665650515e-05,-5.1401253665650515e-05,1.0000000000000000e+00
+60,4.5116286702089390e+01,1.3635660242474829e-01,3.8474842126758285e-01,5.0941896089957439e-01,1.0076775625847574e-01,-6.6588419265878676e-03,-3.2292507987109829e-02,9.8989518395093545e-01
+99,6.0941402944749235e+01,5.6329632662556939e-01,1.6508974083429478e-01,2.8363900072384202e-01,1.4600652272940820e-01,-9.4503553504148817e-03,-8.1107979806100003e-02,9.5744613494697117e-01
+149,3.1250050187602852e+01,9.2019661964820609e-01,3.7783373861371913e-03,9.2371464543242410e-02,7.3204243866698743e-02,-4.9622131664648879e-04,-4.4722231624103201e-02,9.3101182228834278e-01
+199,1.3831174797835914e+01,9.6871352412464840e-01,6.5656925040025985e-04,3.5975407531667690e-02,3.0054736574624694e-02,3.9747232336301581e-03,1.9073206439169550e-02,9.4072464272716960e-01
+299,2.4009864451104130e+00,9.9499349091655198e-01,9.9584290598387748e-05,1.7437058100832695e-02,4.8139270752370631e-03,4.2277664579140582e-03,5.0009060755114958e-02,9.8950899148243476e-01
+349,9.4836453242079399e-01,9.9796063893219689e-01,5.2352355165063066e-05,1.6198462451437204e-02,1.9377442578147675e-03,1.7892857982496535e-04,1.9617772812423967e-02,9.9999574382061696e-01
+449,1.5914313124945889e-01,9.9961814698125151e-01,2.5461006769429786e-05,1.5334951073125025e-02,3.3144008551800805e-04,4.3808552516745440e-27,3.2616629597560521e-04,1.0000000000000000e+00
+599,1.6805285673456662e-02,9.9992369564595407e-01,2.0663311129012729e-05,1.5203845072689088e-02,3.5010987338793960e-05,5.3073494937017952e-61,6.9923499097932138e-07,1.0000000000000000e+00
diff --git a/tests/unitTests/reference_generators/active_stress/land_niederer/README.md b/tests/unitTests/reference_generators/active_stress/land_niederer/README.md
new file mode 100644
index 000000000..32939fa77
--- /dev/null
+++ b/tests/unitTests/reference_generators/active_stress/land_niederer/README.md
@@ -0,0 +1,24 @@
+# Land-Niederer active-stress reference
+
+The committed reference trajectory was produced with KB Ko's [`generate_ref_solution`](https://github.com/kko27/generate_ref_solution.git) using the Land-Niederer 
+human ventricular cardiomyocyte parameter set for multiscale simulations. It 
+was not generated by svMultiPhysics. 
+
+The reference uses 600 outer updates with `dt=1 ms`. At the start of each
+update, calcium follows the normalized double-exponential transient prescribed
+by `ActiveStressTrajectoryTest`: baseline `1e-4 mM`, peak `9e-4 mM`, rise and
+decay constants `20 ms` and `50 ms`, and onset `10 ms`. Sarcomere length follows
+the same raised-cosine trajectory, shortening from `2.2 um` at `30 ms` to
+`2.134 um` at `150 ms` and recovering to `2.2 um` at `350 ms`. Its rate is the
+forward difference `[SL(t+dt)-SL(t)]/dt`.
+
+For each update, the authors' `solve_time_step` receives time in seconds,
+calcium in micromolar, sarcomere length in micrometers, and its rate in
+micrometers per second. The authors' code advances the states with Forward-Euler 
+time steps. The Active tension is divided by 1000 to convert kPa to MPa.
+
+The committed CSV stores `step,Ta,s0,...,s6` at update indices 
+`0,10,30,60,99,149,199,299,349,449,599`, where index `N` is the result of the
+update starting at `t=N ms`. No ready-to-run generator is included because 
+reproducing the data depends on installation of MATLAB. The trusted data remain in
+`tests/unitTests/reference_data/active_stress_land_niederer_twitch.csv`.