Add RDQ20-MF active stress model

[samibismar]

Problem

svMultiPhysics currently does not include the RDQ20-MF mean-field active mechanics model described by Regazzoni, Dedè, and Quarteroni (2020).

Adding this model would allow simulations using a biophysically detailed sarcomere-level active stress formulation that captures calcium activation, length-dependent activation, and force–velocity behavior.

Solution

Add RDQ20-MF as a new active stress model within the ActiveStress framework introduced in PR #578.

The implementation should include:

• the 20-state RDQ20-MF mean-field dynamics,
• model-specific parameter handling through the XML system,
• integration with the active stress model factory,
• MPI parameter distribution,
• active tension calculation,
• and an electromechanics regression test case.

The model formulation is based on:

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

Additional context

The implementation was developed using the authors' reference implementation as a comparison source:

https://github.com/FrancescoRegazzoni/cardiac-activation

Validation was performed by comparing RDQ20-MF state evolution and active tension against the reference implementation across multiple calcium, sarcomere length, and velocity protocols.

Code of Conduct

• I agree to follow this project's Code of Conduct and Contributing Guidelines

[samibismar]

Implemented the RDQ20-MF active stress model on the following branch:

The implementation includes:

• a new RDQ20-MF active stress model following the ActiveStress framework,
• the 20-state mean-field dynamics (16 regulatory-unit states + 4 crossbridge moments),
• XML parameter integration and MPI parameter distribution,
• active tension computation,
• and a slab electromechanics regression test case.

The model was implemented based on the formulation described in Regazzoni et al. (2020) and compared against the authors' C++ reference implementation.

Validation performed so far:

• RDQ20-MF state evolution and active tension compared against the reference implementation across calcium, sarcomere length, and velocity protocols.
• Long-horizon comparisons of all 20 states showed agreement at machine precision.
• Additional force–calcium, length-dependent activation, and force–velocity checks reproduce the expected model behavior.

One additional mechanics question came up during validation regarding how the scalar active tension is assembled into the continuum stress tensor. I discussed this with @michelebucelli, and this appears to be more closely related to the mechanics assembly layer than the RDQ20-MF model itself.

The paper formulation uses the normalized (1/lambda_f) factor in the active stress expression, while the current svMultiPhysics mechanics assembly through compute_pk2cc follows the existing second Piola-Kirchhoff stress pathway.

I would appreciate discussion on whether this difference reflects an intentional mechanics convention in svMultiPhysics, or whether the active tension assembly should be adjusted to account for the selected stress measure.

[michelebucelli]

I discussed this with @michelebucelli, and this appears to be more closely related to the mechanics assembly layer than the RDQ20-MF model itself. The paper formulation uses the normalized (1/lambda_f) factor in the active stress expression, while the current svMultiPhysics mechanics assembly through compute_pk2cc follows the existing second Piola-Kirchhoff stress pathway.

More in detail: the code currently implements active stress by adding the following to the second Piola-Kirchhoff stress tensor:

$$ S_\text{act} = T_\text{act} \mathbf{f}_0 \otimes \mathbf{f}_0;. $$

This is inherited from the old implementation of fiber reinforcement stress (but honestly I hadn't given too much thought to this when working on ActiveStress classes).

What is often found in publications, instead, is the following term added to the first Piola-Kirchhoff stress tensor:

$$ P_\text{act} = T_\text{act} \frac{F\mathbf{f}_0 \otimes \mathbf{f}_0}{|F\mathbf{f}_0|};. $$

The two formulations are not equivalent. Given the relation $P_\text{act} = F S_\text{act}$, they differ by a factor $1 / |F\mathbf{f}_0| = 1 / \lambda_f$, which is what @samibismar noticed.

I lean towards the second formulation, because it seems to me that the active tension has a better physical interpretation in that case (force in the deformed configuration, instead of force in the reference configuration). Still, I haven't (yet) surveyed the literature of cardiac electromechanics about this, so I might be wrong. I'd also like to hear what other people think @aabrown100-git @dseyler @kko27.

[kko27]

@michelebucelli I think what you are asking is that we have two choices to define the active contribution to the 2nd PK stress $S_\text{act}$:

1. $S_\text{act} = T_\text{act} \mathbf{f}_0 \otimes \mathbf{f}_0$ or
2. $S_\text{act} = \frac{T_\text{act}}{\lambda_f} \mathbf{f}_0 \otimes \mathbf{f}_0$

and we need to decide between these two. If we adopt the second option, the deformation dependence of $1/\lambda_f$ needs to be included in the consistent tangent, as noted in @vvedula22's 3rd point in issue #523.

I've started to survey the literature and there doesn't seem to be a convention. Wong (2013) using Nash-Panfilov activation and Myklebust (2025) using Land-Niederer activation do not normalize (Option 1). Whereas Levrero-Florencio (2020) and Franzone (2018) normalize (though slightly differently to the options listed above).

If you look at the mechanics benchmarks papers, I think the non-normalized option is used: 2025 benchmark and 2015 benchmark. But these benchmark papers use a prescribed active stress waveform. This may explain why we have currently have option 1 in our codebase. To my knowledge, there is not yet a coupled cardiac electromechanics benchmark with cellular active-tension models. So I'm not sure if a consensus exists yet, especially for length-dependent active stress models. I'm wondering if this is then a cellular model-specific choice...

Tagging @javijv4 to get additional thoughts

[aabrown100-git]

Very interesting discussion, I'm not too sure which option is better. One thought though, it seems like this choice would influence the degree of length-dependent activation. So if the RDQ20-MF model accurately models length-dependent activation and uses option 2, I think option 2 is a better choice.

If both options are seen in the literature, perhaps we can support both and choose depending on the specific active stress model. I don't think the tangent contribution for option 2 will be that hard to include?

[michelebucelli]

Thank you both for your feedback (and sorry for missing the comment about this in #523)!

I have compiled the references posted by @kko27 and some others into the table reported below (in the collapsible section; feel free to edit this comment to add more if you want).

Table with different normalization strategies in different papers

paper	DOI	notes	T	T/lambda	T/lambda2	JT/lambda2
Kerchoffs et al. (2003)	https://doi.org/10.1114/1.1566447	see eq. (5)	x
Wong et al. (2013)	https://doi.org/10.1002/cnm.2565	see eq. (27); also applies to other papers from the Kuhl group	x
Land et al. (2015)	https://doi.org/10.1098/rspa.2015.0641	see eq. (2.13)	x
Myklebust et al. (2025)	https://doi.org/10.1016/j.cma.2025.118270	see eq. (6) and (20)	x
Arostica et al. (2025)	https://doi.org/10.1016/j.cma.2024.117485	see eq. (2) and (5)	x
Santiago et al. (2018)	https://doi.org/10.1002/cnm.3140	see eq. (6) and (13); cites Humphrey (2003) as source (I don't have access to that one)	x
Ambrosi, Pezzuto (2012)	https://doi.org/10.1007/s10659-011-9351-4	see eq. (4.10); the paper also reviews other choices, but calls the normalized one a "physiologically motivated" one (section 4.2); this is possibly the original source for the iHEART choice		x
Fedele et al. (2023)	https://doi.org/10.1016/j.cma.2023.115983	see eq. (15a); also applies to other papers from the iHEART group		x
Gerach et al. (2021)	https://doi.org/10.3390/math9111247	see eq. (6)		x
Levrero-Florencio et al. (2020)	https://doi.org/10.1016/j.cma.2019.112762	see eq. (2.11)			x
Karabelas et al. (2018)	https://doi.org/10.3389/fphys.2018.00538	see eq. (22); they cite Pathmanathan, Whiteley (2009) as original source			x
Pathmanathan, Whiteley (2009)	https://doi.org/10.1007/s10439-009-9663-8	see eq. (8)			x
Augustin et al. (2019)	https://doi.org/10.1007/s10237-019-01268-5	see eq. (9)			x
Colli Franzone et al. (2018)	https://doi.org/10.3389/fphys.2018.00268	see unnumbered eq. after eq. (4)				x
Land et al. (2017)	https://doi.org/10.1016/j.yjmcc.2017.03.008	unspecified; possibly inherited by Strocchi et al. (2023)	?	?	?	?
Strocchi et al. (2023)	https://doi.org/10.1371/journal.pcbi.1011257	unspecified; points to Land et al. (2017) as a source	?	?	?	?

As you point out, it seems that there is no clear consensus about the normalization strategy: some have no normalization at all, some normalize by the fiber stretch, some by the square of the fiber stretch.

My instinct (but I haven't thought deeply enough about this, and my continuum mechanics is a bit rusty) is that some normalization better, based on the following understanding (which I need to double-check):

1. without normalization, $T_\text{act}$ is the stress per unit reference area, along reference fibers $\mathbf{f}_0$;
2. normalizing by $\lambda$, $T_\text{act}$ is the stress per unit deformed area, along reference fibers $\mathbf{f}_0$;
3. normalizing by $\lambda^2$, $T_\text{act}$ is the stress per unit deformed area, directed along current fibers $F\mathbf{f}_0$.

In my mind, interpretation 3 is better than 2 is better than 1.

Regardless, since all of these are present in the literature, I think it makes sense to have a parameter to switch between the three different strategies (which would also facilitate understanding how big of a difference they make).

[samibismar]

2. normalizing by $\lambda$, $T_\text{act}$ is the stress per unit deformed area, along reference fibers $\mathbf{f}_0$

For the $\frac{T_{act}}{\lambda_f}$ case, would $T_{act}$ instead be interpreted as a force directed along the current unit fiber $f=\frac{Ff_0}{\lambda_f}$, per unit reference area, since

$$ P_{act}=FS_{act}=T_{act}f \otimes f_0 ?$$

[michelebucelli]

@samibismar Another way of seeing my argument is the following (denoting by $S_\text{act}$ the second PK tensor, by $P_\text{act}$ the first PK tensor, and by $\sigma_\text{act}$ the Cauchy tensor).

Without normalization,

$$ \left . \begin{aligned} & S_\text{act} = T_\text{act} \mathbf{f}_0 \otimes \mathbf{f}_0 \\ & \Vert \mathbf{f}_0 \otimes \mathbf{f}_0 \Vert = 1 \end{aligned} \right\rbrace \implies \Vert S_\text{act} \Vert = T_\text{act} $$

Normalizing by $1/\lambda_f$,

$$ \left . \begin{aligned} & S_\text{act} = T_\text{act} \frac{\mathbf{f}_0 \otimes \mathbf{f}_0}{\lambda_f} \\ & P_\text{act} = F S_\text{act} = T_\text{act} \frac{F\mathbf{f}_0 \otimes \mathbf{f}_0}{\lambda_f} \\ & \left\Vert \frac{F\mathbf{f}_0 \otimes \mathbf{f}_0}{\lambda_f} \right\Vert = 1 \end{aligned} \right\rbrace \implies \Vert P_\text{act} \Vert = T_\text{act} $$

Normalizing by $1/\lambda_f^2$,

$$ \left . \begin{aligned} & S_\text{act} = T_\text{act} \frac{\mathbf{f}_0 \otimes \mathbf{f}_0}{\lambda_f^2} \\ & \sigma_\text{act} = \frac{1}{J} F S_\text{act} F^T = \frac{1}{J} T_\text{act} \frac{F\mathbf{f}_0 \otimes F\mathbf{f}_0}{\lambda_f^2} \\ & \left\Vert \frac{F\mathbf{f}_0 \otimes F\mathbf{f}_0}{\lambda_f^2} \right\Vert = 1 \end{aligned} \right\rbrace \implies \Vert \sigma_\text{act} \Vert = \frac{T_\text{act}}{J} $$

Normalizing by $J/\lambda_f^2$, like in Colli Franzone et al. (2018),

$$ \left . \begin{aligned} & S_\text{act} = J T_\text{act} \frac{\mathbf{f}_0 \otimes \mathbf{f}_0}{\lambda_f^2} \\ & \sigma_\text{act} = \frac{1}{J} F S_\text{act} F^T = T_\text{act} \frac{F\mathbf{f}_0 \otimes F\mathbf{f}_0}{\lambda_f^2} \\ & \left\Vert \frac{F\mathbf{f}_0 \otimes F\mathbf{f}_0}{\lambda_f^2} \right\Vert = 1 \end{aligned} \right\rbrace \implies \Vert \sigma_\text{act} \Vert = T_\text{act} $$

If all of this is correct, then the different choices correspond to interpreting $T_\text{act}$ as the norm of a different stress tensor. Accordingly, their physical interpretations differ in a way that mirrors the interpretations of the stress tensors, which is what my previous comment was gesturing at (perhaps not very clearly).

Now that I write it like this, options 1, 2 and 4 seem different but legitimate. Option 3 seems a more unusual hybrid ($T_\text{act}$ is not the norm of any of the three stress tensors), but perhaps it still makes sense in the incompressible limit (where $J \approx 1$).

[javijv4]

Hi everyone,

I looked at some continuum mechanics references, and I believe the current implementation and option 1 of the last @michelebucelli comment come from an energy formulation, which is how most constitutive relationships are developed.

If we describe the energy of a material as

$$ \Psi(\mathbf{C}) = \Psi_{pas}(\mathbf{C}) + \Psi_{act}(I_{4f}) $$

where $\Psi_{pas}$ and $\Psi_{act}$ are the passive and active contributions. Note that I assumed the active term only depends on the deformation across the fiber, which is usually the case.

Then, the active PK2 stress $\mathbf{S}_{act}$ can be derived as

$$ \mathbf{S}_{act} = 2\frac{\partial \Psi_{act}}{\partial \mathbf{C}} = 2\frac{\partial \Psi_{act}}{\partial I_{4f}}\frac{\partial I_{4f}}{\partial \mathbf{C}} = 2\Psi_{act}'(I_{4f})\mathbf{f}_0 \otimes \mathbf{f}_0 $$

using $\partial I_{4f}/\partial \mathbf{C} = \mathbf{f}_0\otimes\mathbf{f}_0$.

From where we get the

$$ \boldsymbol{\sigma}_{act} = \frac{2\Psi_{act}'(I_{4f})}{J} \mathbf{f} \otimes \mathbf{f} $$

And we can replace $2\Psi_{act}'(I_{4f})$ with $\sigma_a$ to avoid defining the energy term.

For that reason, I am leaning toward just leaving the current implementation. However, since other papers have used other variations adding the ability to choose would also be okay with me, but if we do that, we should be careful not to confuse users.

[michelebucelli]

Thank you @javijv4!

Within this energy-based derivation, I suppose one could incorporate the scaling factors into the definition of $\Psi_{act}^\prime$. For example, for the RDQ20-MF model, the code could define

$$T_\text{act} = \frac{\tilde{T}_\text{act}}{\lambda_f}$$

where $\tilde T_{act}$ is the active tension defined in the original paper. This way, the model would behave consistently with the publications, while maintaining the same definition for the active stress tensor $S_\text{act}$ as in the current implementation.

Do you think this would make sense?

[javijv4]

@michelebucelli, yes! I like that idea. Also, that form can be derived from an energy form, and it ensures its convexity (at least when $\tilde T_{act}$ is constant), which is not really the case when we don't normalize by $\lambda_f$. See Ambrosi & Pezzuto's paper. I guess that is a point in favor of always normalizing by $\lambda_f$, but I like the idea of leaving it open to the model and in the context of different definitions of $\Psi_{act}'$, such that we can reproduce the different literature.

[michelebucelli]

Perfect!

@samibismar Then I suggest you proceed by incorporating the normalization into RDQ20MF::compute_active_tension_local. We can revise this decision down the line if necessary.

After that, I looked at the code a bit and it seems generally good, so I think you can open a PR for review. Thanks!