Add comparison with non-frankstartlig isometric twitch

Author: finsberg

demo/time_dependent/frank_starling_twitch.py

@@ -7,26 +7,63 @@
 # 1. Stretch the tissue to a specific length and lock the boundaries (isometric condition).
 # 2. Measure the **passive** stress generated by stretching the hyperelastic collagen/myocardium matrix.
 # 3. Activate the tissue and measure the additional **active** stress generated.
-# 4. Repeat this for different stretch lengths to reconstruct the ascending limb of the cardiac length-tension curve.
-
+# 4. Compare a standard active stress model against our custom Frank-Starling model.
+#
+# ## Mathematical Formulation of Active Stress
+# In a standard active stress approach, the total 2nd Piola-Kirchhoff active stress acts purely along the reference fiber direction $\mathbf{f}_0$ and is independent of stretch:
+#
+# $$
+# \mathbf{S}_{\text{act}} = a(t) \, \mathbf{f}_0 \otimes \mathbf{f}_0
+# $$
+#
+# Where $a(t)$ is the time-dependent activation scalar.
+#
+# For the **Frank-Starling** mechanism, we introduce a stretch-dependent scalar multiplier $g(\lambda)$. The total active First Piola-Kirchhoff stress ($\mathbf{P}_{\text{act}}$) evaluated by the FEM solver becomes:
+#
+# $$
+# \mathbf{P}_{\text{act}} = \mathbf{F} \mathbf{S}_{\text{act}} = \mathbf{F} \Big( \underbrace{a(t) \cdot g(\lambda)}_{T_a} \, \mathbf{f}_0 \otimes \mathbf{f}_0 \Big)
+# $$
+#
+# Where $\mathbf{F}$ is the deformation gradient, and $\lambda = \sqrt{\mathbf{f}_0 \cdot (\mathbf{F}^T \mathbf{F} \mathbf{f}_0)}$ is the dynamic fiber stretch.
+#
+# The multiplier $g(\lambda)$ is defined as a piecewise function representing the ascending limb of the length-tension curve:
+#
+# $$
+# g(\lambda) =
+# \begin{cases}
+#     a_{\min} & \text{if } \lambda \le \lambda_{\text{thres}} \\
+#     a_{\min} + m (\lambda - \lambda_{\text{thres}}) & \text{if } \lambda_{\text{thres}} < \lambda \le \lambda_{\text{opt}} \\
+#     a_{\max} & \text{if } \lambda > \lambda_{\text{opt}}
+# \end{cases}
+# $$
+#
+# Where the slope $m$ is calculated as:
+#
+# $$
+# m = \frac{a_{\max} - a_{\min}}{\lambda_{\text{opt}} - \lambda_{\text{thres}}}
+# $$
+#
 # %% [markdown]
 # ## 1. Imports and Problem Setup
-# First, we import the necessary modules. We use `pulse` for the mechanics framework, `dolfinx` for the finite element backend, and `mpi4py` for parallel execution.
+# First, we import the necessary modules. We use `pulse` for the mechanics framework, `dolfinx` for the finite element backend, `mpi4py` for parallel execution, and `matplotlib` for plotting our final results.
 
 # %%
 from mpi4py import MPI
 import dolfinx
 import pulse
 import numpy as np
 import ufl
+import matplotlib.pyplot as plt
 from scifem import evaluate_function
 
 # %% [markdown]
 # ## 2. Defining the Isometric Twitch Experiment
 # We define a function that takes a prescribed pre-stretch (in mm) and an active twitch tension (in kPa). It will stretch a 10x1x1 mm block of tissue, lock it in place, and then apply the active tension.
+#
+# We include a boolean flag `use_frank_starling` to easily toggle between the standard active stress model and our stretch-dependent model.
 
 # %%
-def run_isometric_twitch(pre_stretch_mm: float, twitch_activation: float):
+def run_isometric_twitch(pre_stretch_mm: float, twitch_activation: float, use_frank_starling: bool = True):
     """
     Runs an isometric contraction on a 10x1x1 slab.
     1. Stretches the right face by `pre_stretch_mm` and locks it.
@@ -50,14 +87,17 @@ def run_isometric_twitch(pre_stretch_mm: float, twitch_activation: float):
         dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(0.0)), "kPa",
     )
 
-    # Set up the passive model
+    # Set up the passive material model
     material_params = pulse.HolzapfelOgden.transversely_isotropic_parameters()
     passive_model = pulse.HolzapfelOgden(
         f0=f0, s0=s0, **material_params,
     )
 
-    # Set up our custom stretch-dependent active model
-    active_model = pulse.FrankStarlingActiveStress(f0=f0, activation=Ta)
+    # Toggle between the Frank-Starling model and the standard Active Stress model
+    if not use_frank_starling:
+        active_model = pulse.ActiveStress(f0=f0, activation=Ta)
+    else:
+        active_model = pulse.FrankStarlingActiveStress(f0=f0, activation=Ta)
 
     # Combine into a fully compressible cardiac model
     model = pulse.CardiacModel(
@@ -107,13 +147,16 @@ def dirichlet_bc(
     # Set up the FEniCSx-Pulse mechanics problem
     parameters = {"mesh_unit": "mm"}
     problem = pulse.StaticProblem(model=model, geometry=geo, bcs=bcs, parameters=parameters)
+
+    # CRITICAL: Register the displacement field for the Frank-Starling stretch tracking
+    if use_frank_starling:
+        active_model.register(problem.u)
+
     point = np.array([[L, 0.5, 0.5]])
 
     # --- Phase 1: Passive Pre-Stretch ---
     # Solve for the passive stretch equilibrium
     problem.solve()
-    u_pre_stretch = evaluate_function(problem.u, point)[0][0] # X-displacement at the tip
-    length_pre_stretch = L + u_pre_stretch
 
     # Formulate the average fiber stress
     F = ufl.variable(ufl.grad(problem.u) + ufl.Identity(3))
@@ -136,9 +179,6 @@ def dirichlet_bc(
         Ta.assign(twitch_activation * (step + 1) / nsteps)
         problem.solve()
 
-    u_twitch = evaluate_function(problem.u, point)[0][0]
-    length_twitch = L + u_twitch
-
     # Calculate total stress after activation
     total_force = mesh.comm.allreduce(dolfinx.fem.assemble_scalar(Tf), op=MPI.SUM) / volume
 
@@ -147,29 +187,90 @@ def dirichlet_bc(
 
     return passive_force, active_force
 
-
 # %% [markdown]
-# ## 3. Running the Curve
-# To generate the length-tension curve, we loop over increasing amounts of stretch (from 0% up to 15%) and run the experiment. We can then observe the non-linear stiffening of the passive matrix against the dynamic active response of the tissue.
+# ## 3. Running the Curve and Comparison
+# We will loop over increasing amounts of stretch (from 0% up to 15%) and run the experiment twice per stretch: once with the standard active stress, and once with the Frank-Starling model.
 
 # %%
-if __name__ == "__main__":
-    twitch_force = 1.0 # kPa of active baseline tension
 
-    # We will test stretches from 0 mm (0%) up to 1.5 mm (15%)
-    stretch_amounts = [0.0, 0.25, 0.5, 0.75, 1.0, 1.25, 1.5]
+twitch_force = 1.0 # kPa of active baseline tension
+
+# Test stretches from 0 mm (0%) up to 1.5 mm (15%)
+stretch_amounts = [0.0, 0.25, 0.5, 0.75, 1.0, 1.25, 1.5]
+strain_pcts = [(s / 10.0) * 100 for s in stretch_amounts]
+
+passive_stresses = []
+active_stresses_std = []
+active_stresses_fs = []
+
+print(f"{'Stretch (%)':>12} | {'Passive Stress':>18} | {'Std Active Stress':>20} | {'FS Active Stress':>20}")
+print("-" * 78)
+
+for s, pct in zip(stretch_amounts, strain_pcts):
+    # 1. Run Standard Model
+    p_force, a_force_std = run_isometric_twitch(s, twitch_force, use_frank_starling=False)
+
+    # 2. Run Frank-Starling Model
+    _, a_force_fs = run_isometric_twitch(s, twitch_force, use_frank_starling=True)
 
-    print(f"{'Stretch (mm)':>12} | {'Stretch (%)':>12} | {'Passive Stress (kPa)':>20} | {'Active Stress (kPa)':>20}")
-    print("-" * 75)
+    # Store results
+    passive_stresses.append(p_force)
+    active_stresses_std.append(a_force_std)
+    active_stresses_fs.append(a_force_fs)
 
-    for stretch in stretch_amounts:
-        p_force, a_force = run_isometric_twitch(pre_stretch_mm=stretch, twitch_activation=twitch_force)
-        strain_pct = (stretch / 10.0) * 100
-        print(f"{stretch:12.2f} | {strain_pct:12.2f} | {p_force:20.4f} | {a_force:20.4f}")
+    print(f"{pct:12.2f} | {p_force:18.4f} | {a_force_std:20.4f} | {a_force_fs:20.4f}")
+
+# %% [markdown]
+# ## 4. Plotting the Results
+# While a table is useful, plotting the results makes the physiological phenomenon instantly clear. We will use `matplotlib` to plot the Active Stress (comparing both models) on the left axis, and the Passive hyperelastic stress on the right axis.
+
+# %%
+fig, ax1 = plt.subplots(figsize=(9, 6))
+
+# Plot Active Stresses on the primary y-axis
+ax1.set_xlabel('Stretch (%)', fontsize=12)
+ax1.set_ylabel('Active Stress (kPa)', color='tab:blue', fontsize=12)
+
+ax1.plot(
+    strain_pcts, active_stresses_fs, marker='o', linewidth=2, color='tab:blue',
+    label='Frank-Starling Active Stress',
+)
+ax1.plot(
+    strain_pcts, active_stresses_std, marker='s', linestyle='--', linewidth=2, color='tab:cyan',
+    label='Standard Active Stress',
+)
+
+ax1.tick_params(axis='y', labelcolor='tab:blue')
+ax1.grid(True, linestyle='--', alpha=0.6)
+
+# Create a twin axis for the Passive Stress
+ax2 = ax1.twinx()
+ax2.set_ylabel('Passive Stress (kPa)', color='tab:red', fontsize=12)
+ax2.plot(
+    strain_pcts, passive_stresses, marker='^', linewidth=2, color='tab:red',
+    label='Passive Stress (Holzapfel-Ogden)',
+)
+ax2.tick_params(axis='y', labelcolor='tab:red')
+
+# Add legends
+lines_1, labels_1 = ax1.get_legend_handles_labels()
+lines_2, labels_2 = ax2.get_legend_handles_labels()
+ax1.legend(lines_1 + lines_2, labels_1 + labels_2, loc='upper center', bbox_to_anchor=(0.5, -0.15), ncol=3)
+
+plt.title('Isometric Contraction: Frank-Starling vs. Standard Active Stress', fontsize=14)
+fig.tight_layout()
+plt.show()
 
 # %% [markdown]
 # ### Conclusion
 #
-# The output from the script prints a table showcasing:
-# 1. The **Passive Stress** climbs exponentially as the tissue is stretched (capturing the non-linear stiffening of the collagen matrix).
-# 2. The **Active Stress** scales independently and mimics the Frank-Starling law. At 0% stretch, the tissue produces low active force. As it is stretched towards 15%, the local fibers dynamically detect the stretch and generate a much larger active contractile stress.
+# The plot generated by the script perfectly illustrates two primary mechanics of the heart wall:
+# 1. **Exponential Passive Stiffening (Red Line):** As the tissue is stretched past 10%, the collagen fibers un-crimp and become highly resistant to further stretch.
+# 2. **Frank-Starling Law (Dark Blue vs Light Blue):** The standard model (dashed light blue) assumes a constant, 100% force capacity regardless of how compressed or stretched the sarcomeres are. The Frank-Starling model (solid dark blue) dynamically scales the force—starting lower at rest, and sharply climbing as the tissue stretches toward its optimal operating length.
+#
+# **Wait, why does the Standard Active Stress slightly increase?**
+# You might notice that the standard model's active stress isn't perfectly flat across stretches. This is actually physically and mathematically correct! In continuum mechanics, the active stress is defined in the *reference* configuration (2nd Piola-Kirchhoff stress). When we measure the *true* (Cauchy) stress in the *deformed* configuration, the solver applies a geometric "push-forward".
+#
+# As the tissue stretches, it physically narrows due to incompressibility (the Poisson effect). Because the cross-sectional area shrinks, the true stress (Force / Area) mechanically increases, even if the chemical contraction capability remained identical.
+#
+# Thus, the Standard model increases purely due to **kinematic/geometric stiffening**, while the Frank-Starling model increases due to that same geometric effect *plus* the **biological recruitment** of more myosin-actin cross-bridges.

src/pulse/active_stress.py

@@ -332,7 +332,7 @@ def transversely_active_stress_strain_energy(Ta, C, f0, eta=0.0):
 
     Arguments
     ---------
-    Ta : dolfinx.fem.Function or dolfinx.femConstant
+    Ta : dolfinx.fem.Function or dolfinx.fem.Constant
         A scalar function representing the magnitude of the
         active stress in the reference configuration (first Piola)
     C : ufl.Form
@@ -364,7 +364,7 @@ def transversely_active_stress(Ta, f0, eta=0.0):
 
     Arguments
     ---------
-    Ta : dolfinx.fem.Function or dolfinx.femConstant
+    Ta : dolfinx.fem.Function or dolfinx.fem.Constant
         A scalar function representing the magnitude of the
         active stress in the reference configuration (first Piola)
     f0 : dolfinx.fem.Function