Add Regazzoni active stress model

[samibismar]

Current situation

This PR adds the RegazzoniActiveStress model, implementing the RDQ20-MF formulation described by Regazzoni, Dedè, and Quarteroni (2020), using the ActiveStress framework introduced in PR #578.

The implementation includes:

• the 20-state RDQ20-MF system, consisting of 16 regulatory-unit probabilities and 4 crossbridge moments;
• model-specific XML parameter parsing and MPI distribution;
• registration with the ActiveStress model factory;
• calcium-, sarcomere-length-, and velocity-dependent state evolution;
• active-tension computation;
• and an electromechanics regression test.

The local active-stress interface is extended to provide fiber stretch to each model. RDQ20-MF first computes the paper tension

$$ \widetilde T_{\mathrm{act}} = a_{\mathrm{XB}} \left(\mu_P^1+\mu_N^1\right)\phi(SL), $$

then returns

$$ T_{\mathrm{code}} = \frac{\widetilde T_{\mathrm{act}}}{\lambda_f}. $$

This normalization is required because svMultiPhysics assembles the returned coefficient as a second Piola–Kirchhoff stress contribution:

$$ \mathbf S_{\mathrm{act}} = T_{\mathrm{code}}\ \mathbf f_0\otimes\mathbf f_0. $$

After conversion to first Piola–Kirchhoff stress, this gives

$$\mathbf P_{\mathrm{act}} = \widetilde T_{\mathrm{act}}\ \mathbf f\otimes\mathbf f_0, $$

consistent with the RDQ20 paper convention.

Related to #594.

Release Notes

• Add the RDQ20-MF active-stress model.
• Add RDQ20-MF XML parameters and MPI parameter distribution.
• Extend the local active-stress interface to receive fiber stretch.
• Normalize the RDQ20-MF returned coefficient for the existing PK2 mechanics assembly.
• Add the slab_RDQ20MF electromechanics regression case.

Documentation

The RDQ20-MF class documentation describes:

• the model states and parameters
• the distinction between the paper tension $\widetilde T_{\mathrm{act}}$ and the returned PK2 coefficient $T_{\mathrm{code}}$
• the fiber-stretch normalization
• and the exception raised for zero or near-zero fiber stretch

The slab_RDQ20MF README documents:

• the coupled electromechanics setup
• the svMultiPhysics fiber/sheet/sheet-normal directional weighting
• the meaning of the stored Active_tension_* fields
• and the reference-validation procedure

Testing

• The project build completed successfully.
• The slab_RDQ20MF pytest case passed for all three configured process-count variants.
• The RDQ20-MF state evolution and paper tension were validated against the authors' C++ reference implementation at commit 26f05df.
• After normalization, the reconstructed paper tension

$$ \widetilde{T}_{\mathrm{act}} = \lambda_f T_{\mathrm{code}} $$

was compared node-by-node against the same reference implementation and agreed to approximately machine precision.

• Additional checks at $\lambda_f = 0.8$, $1.0$, and $1.2$ confirmed

$$ \lambda_f T_{\mathrm{code}} = \widetilde{T}_{\mathrm{act}} $$

to floating-point precision.

• A local full-suite run completed with 206 tests passing and 27 skipped. The skipped tests require optional PETSc or Trilinos build configurations.
• Twelve sv0D coupling tests could not run because the local svZeroDSolver shared library was unavailable. The same twelve tests fail with the same missing-library error on main in the same environment.

References

Regazzoni, F., Dede', L., & Quarteroni, A.
“Biophysically detailed mathematical models of multiscale cardiac active mechanics.”
PLOS Computational Biology, 2020.
https://doi.org/10.1371/journal.pcbi.1008294

Code of Conduct & Contributing Guidelines

• I agree to follow the Code of Conduct and Contributing Guidelines.

[claude]

Claude Code Review

This pull request is from a fork — automated review is disabled. A repository maintainer can comment @claude review to run a one-time review.

[codecov]

Codecov Report

❌ Patch coverage is 96.80365% with 7 lines in your changes missing coverage. Please review.
✅ Project coverage is 72.72%. Comparing base (b5fbed0) to head (cb56643).
⚠️ Report is 1 commits behind head on main.

Files with missing lines	Patch %	Lines
Code/Source/solver/active_stress_regazzoni.cpp	97.34%	5 Missing ⚠️
Code/Source/solver/active_stress_regazzoni.h	96.42%	1 Missing ⚠️
Code/Source/solver/active_stress_uniform_steady.h	0.00%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #596      +/-   ##
==========================================
+ Coverage   72.54%   72.72%   +0.17%     
==========================================
  Files         252      255       +3     
  Lines       39032    39285     +253     
  Branches     6678     6731      +53     
==========================================
+ Hits        28317    28571     +254     
+ Misses      10480    10472       -8     
- Partials      235      242       +7     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[ktbolt]

@samibismar May I suggest using a more meaningful name than RDQ20MF ? Maybe just MeanField perhaps qualified by what it actually represents rather than authors and publication date ?

[samibismar]

@ktbolt Would MeanFieldActiveStress be an appropriate class name, with MeanField as the factory label used in the XML input, for example <Model>MeanField</Model>?

I would retain “RDQ20-MF” in the documentation and citations to identify the specific Regazzoni–Dedè–Quarteroni formulation being implemented.

[ktbolt]

@samibismar I think MeanFieldActiveStress is fine and just document it what the implementation is based on.

[michelebucelli]

@samibismar @ktbolt Sorry to comment on this this late, but: I'm not sure that MeanFieldActiveStress is a descriptive name for the model. Mean-field based models are more a general category than a description of this specific model, and the fact that the model is based on the mean-field formalism is, I think, a relatively secondary aspect of it.

As far as I know, the standard way of referring to these cellular models is by the name of the author (and sometimes year, if the same authors have published multiple models), since (I think) it is often impossible to come up with a unique descriptive name. This is the same as with ionic models (and indeed we have the classes AlievPanfilov, BuenoOrovio and so on).

If you feel the acronym RDQ is a bit too obscure, I wouldn't mind using Regazzoni instead (or RegazzoniActiveStress).

[samibismar]

@michelebucelli That makes sense, thanks. I agree that MeanFieldActiveStress is probably too broad if mean-field models are a general category rather than a unique description of this formulation.

I think RegazzoniActiveStress for the C++ class, with Regazzoni as the model/XML label, would fit the existing naming pattern while remaining clear in the active-stress code. I would continue to identify it as the RDQ20-MF formulation of Regazzoni, Dedè, and Quarteroni (2020) in the documentation and validation notes.

@ktbolt, would that naming approach work for you as well?

[michelebucelli]

Thanks @samibismar! Please see my comments below.

[michelebucelli]

I don't really like that this class inherits from ActiveStress rather than ActiveStressODE, even if it is an active stress model based on solving a system of ODEs. Is it because it implements its own custom time-stepping method, or are there other reasons?

I think we should figure out if and how this can be brought under ActiveStressODE (although I can see that this is not exactly trivial). I can see a few ways to do that, and I don't like any of them 😅

1. Generalize the custom time stepping scheme so that it can be moved to ActiveStressODE. This would involve introducing an implicit-explicit scheme, with an interface not unlike that of IonicModel, where state variables are partitioned into gating variables (updated implicitly with a closed-form solution) and others (updated with whatever timestepping method the user chose, preferably explicitly). I don't really like the overhead that this splitting brings.
2. Derive from ActiveStressODE, but then don't use the time stepping infrastructure it provides, and override it with custom time stepping. This would satisfy the "is-one-of" interpretation of inheritance, but break the interface in a number of other ways (e.g. the function getf would make no sense).
3. Implement better explicit time stepping in ActiveStressODE (say, 4th order Runge-Kutta), and see if it works well enough for this model that we can ditch the custom time stepping scheme. This option has a bit of uncertainty related to finding out if RK4 is enough or not.
4. If ActiveStressODE is only used by NashPanfilov, then maybe we can get rid of ActiveStressODE and have each model implement their custom timestepping as needed. I think whether this makes sense or not hinges on whether we expect to use ActiveStressODE in other cases: @kko27 are you using it for the Land-Niederer model, and if so did you need custom timestepping or was forward Euler enough?

I think among these option 3 (the Runge-Kutta one) would be the best, if it works.

[kko27]

@michelebucelli I am using the ActiveStressODE class for my Land-Niederer implementation (for reference, here is my branch) However, I had some convergence issues with the FE time stepping and had to implement RK4 instead.

I am in favor of option 3 (and then 4, if RK4 is not good enough)

[michelebucelli]

Thanks @kko27! I've asked @samibismar to try out the RK4 solver from your branch on the RDQ model, if it works we'll go with option 3.

[samibismar]

I tested the RK4 approach on the Regazzoni model using both an isolated 20-state comparison and the fully coupled electromechanics slab.

In the isolated comparison, stable RK4 timesteps reproduced the refined current-scheme solution under prescribed calcium and fixed kinematics, which suggests that the unified 20-state RHS is consistent with the current model for those imposed inputs.

In the coupled slab:

dt = 1.0 ms: RK4 failed around 11 ms
dt = 0.5 ms: RK4 failed around 11.5 ms
dt = 0.25 ms: RK4 remained stable for the full simulation

The current Regazzoni implementation remained stable at dt = 1.0 ms. This is consistent with the split scheme resolving the faster RU dynamics with smaller internal steps while advancing the XB moments implicitly.

Both unstable coupled runs failed during the early calcium upstroke, with severely unphysical active tension followed by non-finite mechanics fields. The internal RU and XB states were not exported, so I could not determine which state became invalid first.

The result I found is that the shared RK4 approach works for this model only after substantially reducing the outer timestep. I’ve attached a summary figure. The experimental Regazzoni RK4 implementation is available here.