Merge pull request #172 from finsberg/frank-starling

Add frank-starling demo and implementation

Author: finsberg

_toc.yml

@@ -32,6 +32,7 @@ parts:
       - file: "demo/time_dependent/time_dependent_land_circ_lv"
       - file: "demo/time_dependent/time_dependent_land_circ_biv"
       - file: "demo/time_dependent/complete_cycle"
+      - file: "demo/time_dependent/frank_starling_twitch"
     - file: "demo/postprocess_cycle"
 
   - caption: Community

conf.py

@@ -36,6 +36,7 @@
     "sphinx_codeautolink",
     "sphinx_multitoc_numbering",
 ]
+autodoc_typehints = "description"
 external_toc_exclude_missing = True
 external_toc_path = "_toc.yml"
 html_baseurl = ""
@@ -102,3 +103,10 @@
 suppress_warnings = ["mystnb.unknown_mime_type", "bibtex.duplicate_citation"]
 use_jupyterbook_latex = True
 use_multitoc_numbering = True
+nitpick_ignore = [
+    ("py:class", "T"),
+    ("py:class", "pulse.problem.DynamicProblem.Function.T"),
+    ("py:class", "pulse.problem.StaticProblem.Function.T"),
+    ("py:class", "dolfinx.fem.function.Function"),
+    ("py:class", "dolfinx.fem.petsc.NonlinearProblem"),
+]

demo/time_dependent/frank_starling_twitch.py

@@ -0,0 +1,175 @@
+# %% [markdown]
+# # Isometric Contraction and the Frank-Starling Mechanism
+#
+# In cardiac mechanics, the **Frank-Starling mechanism** describes how the heart inherently generates a greater active contraction force when it is subjected to a larger initial stretch (preload). At the cellular level, stretching the tissue moves the sarcomeres into a more optimal overlap configuration, allowing more myosin-actin cross-bridges to form.
+#
+# In this tutorial, we will demonstrate this physiological phenomenon by conducting an **isometric twitch experiment** on a slab of myocardial tissue. We will:
+# 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.
+
+# %% [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.
+
+# %%
+from mpi4py import MPI
+import dolfinx
+import pulse
+import numpy as np
+import ufl
+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.
+
+# %%
+def run_isometric_twitch(pre_stretch_mm: float, twitch_activation: float):
+    """
+    Runs an isometric contraction on a 10x1x1 slab.
+    1. Stretches the right face by `pre_stretch_mm` and locks it.
+    2. Applies `twitch_activation` (kPa) to the active model.
+    3. Returns the average passive and active fiber stresses.
+    """
+    # Create a 10 x 1 x 1 mm slab mesh
+    L = 10.0
+    mesh = dolfinx.mesh.create_box(
+        MPI.COMM_WORLD,
+        [[0.0, 0.0, 0.0], [L, 1.0, 1.0]],
+        [10, 2, 2],
+    )
+
+    # Define microstructure (fibers aligned with the X-axis)
+    f0 = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type((1.0, 0.0, 0.0)))
+    s0 = dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type((0.0, 1.0, 0.0)))
+
+    # Trackers for active tension
+    Ta = pulse.Variable(
+        dolfinx.fem.Constant(mesh, dolfinx.default_scalar_type(0.0)), "kPa",
+    )
+
+    # Set up the passive 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)
+
+    # Combine into a fully compressible cardiac model
+    model = pulse.CardiacModel(
+        material=passive_model,
+        active=active_model,
+        compressibility=pulse.Compressible2(),
+    )
+
+    # --- Geometry and Boundary Conditions ---
+    boundaries = [
+        pulse.Marker(name="X0", marker=1, dim=2, locator=lambda x: np.isclose(x[0], 0)),
+        pulse.Marker(name="X1", marker=2, dim=2, locator=lambda x: np.isclose(x[0], L)),
+    ]
+
+    geo = pulse.Geometry(
+        mesh=mesh,
+        boundaries=boundaries,
+        metadata={"quadrature_degree": 4},
+    )
+
+    def dirichlet_bc(
+        V: dolfinx.fem.FunctionSpace,
+    ) -> list[dolfinx.fem.bcs.DirichletBC]:
+        mesh.topology.create_connectivity(mesh.topology.dim - 1, mesh.topology.dim)
+
+        # 1. Lock the left face
+        facets_fixed = geo.facet_tags.find(1)
+        dofs = dolfinx.fem.locate_dofs_topological(V, 2, facets_fixed)
+        u_fixed = dolfinx.fem.Function(V)
+        u_fixed.x.array[:] = 0.0
+
+        # 2. Stretch the right face in the X direction only
+        facets_stretch = geo.facet_tags.find(2)
+        V_x, _ = V.sub(0).collapse()
+        dofs_x = dolfinx.fem.locate_dofs_topological((V.sub(0), V_x), 2, facets_stretch)
+
+        u_stretch_x = dolfinx.fem.Function(V_x)
+        u_stretch_x.x.array[:] = pre_stretch_mm
+
+        return [
+            dolfinx.fem.dirichletbc(u_stretch_x, dofs_x, V.sub(0)),
+            dolfinx.fem.dirichletbc(u_fixed, dofs),
+        ]
+
+    bcs = pulse.BoundaryConditions(dirichlet=(dirichlet_bc,))
+
+    # Set up the FEniCSx-Pulse mechanics problem
+    parameters = {"mesh_unit": "mm"}
+    problem = pulse.StaticProblem(model=model, geometry=geo, bcs=bcs, parameters=parameters)
+    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))
+    f = F * f0
+    f_norm = f / ufl.sqrt(ufl.inner(f, f))
+
+    # Integrate to find the domain volume
+    volume = mesh.comm.allreduce(
+        dolfinx.fem.assemble_scalar(dolfinx.fem.form(ufl.det(F) * geo.dx)), op=MPI.SUM,
+    )
+
+    # Integrate Cauchy stress projected onto the deformed fiber direction
+    Tf = dolfinx.fem.form(ufl.inner(model.sigma(F) * f_norm, f_norm) * geo.dx)
+    passive_force = mesh.comm.allreduce(dolfinx.fem.assemble_scalar(Tf), op=MPI.SUM) / volume
+
+    # --- Phase 2: Active Twitch ---
+    nsteps = 10
+    for step in range(nsteps):
+        # Ramp up activation over nsteps steps to aid Newton convergence
+        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
+
+    # Isolate the active contribution
+    active_force = total_force - passive_force
+
+    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.
+
+# %%
+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]
+
+    print(f"{'Stretch (mm)':>12} | {'Stretch (%)':>12} | {'Passive Stress (kPa)':>20} | {'Active Stress (kPa)':>20}")
+    print("-" * 75)
+
+    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}")
+
+# %% [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.

docs/api.rst

@@ -76,7 +76,7 @@ geometry
 
 
 problem
-----------------
+-------
 .. automodule:: pulse.problem
     :members:
     :inherited-members:

src/pulse/__init__.py

@@ -29,7 +29,7 @@
     viscoelasticity,
 )
 from .active_model import ActiveModel, Passive
-from .active_stress import ActiveStress
+from .active_stress import ActiveStress, FrankStarlingActiveStress
 from .boundary_conditions import BoundaryConditions, NeumannBC, RobinBC
 from .cardiac_model import CardiacModel
 from .compressibility import (
@@ -106,4 +106,5 @@
     "FixedPointUnloader",
     "PrestressProblem",
     "TargetPressure",
+    "FrankStarlingActiveStress",
 ]

src/pulse/active_model.py

@@ -154,6 +154,9 @@ def P(self, F: ufl.core.expr.Expr, dev: bool = False) -> ufl.core.expr.Expr:
             Cdev = C
         return ufl.diff(self.strain_energy(Cdev), F)
 
+    def register(self, u: dolfinx.fem.Function) -> None:
+        pass
+
 
 class Passive(ActiveModel):
     """Active model with no active component.