Add pseudo-ecg and add niederer benchmark demo

Makefile

@@ -78,6 +78,7 @@ docs: ## generate Sphinx HTML documentation, including API docs
 	jupytext demos/extract_currents.py -o docs/extract_currents.md
 	jupytext demos/tracking_values.py -o docs/tracking_values.md
 	jupytext demos/custom_stimulus_domain.py -o docs/custom_stimulus_domain.md
+	jupytext demos/niederer_benchmark.py -o docs/niederer_benchmark.md
 	jupyter book build -W docs
 
 show-docs:

demos/niederer_benchmark.py

@@ -0,0 +1,178 @@
+# # Niederer benchmark
+#
+# In this example we will use the same setup as in the Niederer benchmark
+# > Niederer SA, Kerfoot E, Benson AP, Bernabeu MO, Bernus O, Bradley C, Cherry EM, Clayton R, Fenton FH, Garny A, Heidenreich E, Land S, Maleckar M, Pathmanathan P, Plank G, Rodríguez JF, Roy I, Sachse FB, Seemann G, Skavhaug O, Smith NP. Verification of cardiac tissue electrophysiology simulators using an N-version benchmark. Philos Trans A Math Phys Eng Sci. 2011 Nov 13;369(1954):4331-51. doi: 10.1098/rsta.2011.0139. PMID: 21969679; PMCID: PMC3263775.
+#
+# However, we will be using an electromechanics model and also compute some pseudo ecg.
+#
+# First me make the necessary imports
+#
+
+from pathlib import Path
+import dolfin
+import matplotlib.pyplot as plt
+import simcardems
+import numpy as np
+
+# And we specify the path to the output directory where we want to store the results
+
+here = Path(__file__).absolute().parent
+outdir = here / "results_niederer_benchmark"
+
+
+# Next we define the stimulus domain
+
+
+def stimulus_domain(mesh):
+    marker = 1
+    subdomain = dolfin.CompiledSubDomain(
+        "x[0] <= L + DOLFIN_EPS && x[1] <= L + DOLFIN_EPS && x[2] <= L + DOLFIN_EPS",
+        L=1.5,
+    )
+    domain = dolfin.MeshFunction("size_t", mesh, mesh.topology().dim())
+    domain.set_all(0)
+    subdomain.mark(domain, marker)
+    return simcardems.geometry.StimulusDomain(domain=domain, marker=marker)
+
+
+# and create the geometry
+
+geo = simcardems.slabgeometry.SlabGeometry(
+    parameters={"lx": 20.0, "ly": 7.0, "lz": 3.0, "dx": 2.0},
+    stimulus_domain=stimulus_domain,
+)
+
+
+# We now create the configurations
+
+config = simcardems.Config(
+    outdir=outdir,
+    coupling_type="fully_coupled_ORdmm_Land",
+    T=1000,
+)
+
+# And create the coupling. Note that, here we are using a different method than usual for creating the coupling, since we need to also supply the geometry.
+
+coupling = simcardems.models.em_model.setup_EM_model_from_config(
+    config=config,
+    geometry=geo,
+)
+
+# Next we create the runner, and solve the problem
+
+runner = simcardems.Runner.from_models(config=config, coupling=coupling)
+runner.solve(T=config.T, save_freq=config.save_freq, show_progress_bar=True)
+
+# Now let us only extract the results of the membrane potential, and compare the values in different points in the mesh. For the slab geometry we can for example evaluate the functions a the minimum and maximum $x$ values and at the center. To do this, we need to first load the results
+
+loader = simcardems.DataLoader(outdir / "results.h5")
+
+# and save the voltage and displacement to xdmf to be viewed in Paraview
+
+simcardems.postprocess.make_xdmffiles(outdir.joinpath("results.h5"), names=["u", "V"])
+
+# <video controls src="./_static/niederer.mp4"></video>
+#
+# Now let us first select some points outside the domain which will serve as our leads for the ecg computation
+
+time_stamps = loader.time_stamps or []
+
+# Here we plot the mesh with the leads as red dots and the ground lead as a black dot
+
+# +
+points = np.array(
+    [
+        (22.0, 8.0, 4.0),
+        (22.0, 0.0, 4.0),
+        (22.0, 8.0, 0.0),
+        (22.0, -1.0, -1.0),
+        (22.0, 3.5, 4.0),
+        (22.0, 3.5, -1.0),
+    ],
+)
+ground = np.array((0.0, -1.0, -1.0))
+fig = plt.figure()
+ax = fig.add_subplot(projection="3d")
+
+dolfin.common.plotting.mplot_mesh(ax, geo.mesh, alpha=0.3)
+ax.scatter(*points.T, color="red")
+ax.scatter(*ground, color="k")
+fig.savefig(outdir / "mesh_with_leads.png")
+
+
+# -
+
+
+# ```{figure} figures/mesh_with_leads.png
+# ---
+# name: mesh_with_leads
+# ---
+# Mesh with leads as red dots and ground lead as black dot
+# ```
+#
+# We will now compute the pseudo-ecg using the recovery formula. In this approach we estimate the extracellular potential using the formula
+#
+# $$
+# \phi_b(\mathbf{x}, t) = \frac{1}{4 \pi \sigma_b} \int_{\Omega} \frac{ M_i \nabla V(\tilde{\mathbf{x}}, t)  (\tilde{\mathbf{x}} - \mathbf{x})}{\| \tilde{\mathbf{x}} - \mathbf{x} \|^3} \mathrm{d}\tilde{\mathbf{x}}
+# $$
+#
+# And we compute the pseudo-ECG by subtracting $\phi_b(\mathbf{x}_p, t) - \phi_b(\mathbf{x}_g, t)$ where $\mathbf{x}_p$ is the position of the lead and $\mathbf{x}_g$ is the position of a ground lead.
+
+
+# +
+def voltage_traces():
+    for t in time_stamps:
+        yield loader.get("ep", "V", t)
+
+
+fig, ax = plt.subplots()
+for point in points:
+    ecg = simcardems.postprocess.ecg(
+        voltage_traces(),
+        sigma_b=1.0,
+        point1=point,
+        point2=(-10.0, -10.0, -10.0),
+    )
+
+    ax.plot(np.array(time_stamps, dtype=float), ecg, label=point)
+ax.legend()
+fig.savefig(outdir / "ecg.png")
+# -
+
+# ```{figure} figures/ecg.png
+# ---
+# name: ecg
+# ---
+# Pseudo-ECG
+# ```
+#
+# Next we use a function for computing the activation map, and indicate that the threshold value should be 0.0. The activation map, will then be a function over the mesh where the value at a given point will be the time it takes for the voltage to first reach 0.0 mV.
+
+activation_map = simcardems.postprocess.activation_map(
+    voltage_traces(),
+    time_stamps=time_stamps,
+    V=loader._functions["ep"]["V"].function_space(),
+    threshold=0.0,
+)
+
+
+# and compute the activation times for all points through the center of the mesh, from $x=0$ to $x=L_x$.
+#
+
+x = np.linspace(0, geo.parameters["lx"], 50)
+act = [
+    activation_map(xi, geo.parameters["ly"] / 2, geo.parameters["lz"] / 2) for xi in x
+]
+fig, ax = plt.subplots()
+ax.plot(x, act)
+ax.set_xlabel("$x$ coordinate")
+ax.set_ylabel("Activation time [ms]")
+ax.grid()
+fig.savefig(outdir / "activation_times.png")
+
+# ```{figure} figures/activation_times_niederer.png
+# ---
+# name: activation_times_niederer
+# ---
+# Activation times through the center of the mesh.
+# ```

docs/_toc.yml

@@ -24,6 +24,7 @@ parts:
     - file: "extract_currents"
     - file: "tracking_values"
     - file: "custom_stimulus_domain"
+    - file: "niederer_benchmark"
     - file: "cli"
     - file: "docker"
     - file: "gui"

src/simcardems/postprocess.py

@@ -14,6 +14,11 @@
 import numpy as np
 import tqdm
 
+try:
+    import ufl_legacy as ufl
+except ImportError:
+    import ufl
+
 from . import utils
 from .datacollector import DataCollector
 from .datacollector import DataGroups
@@ -277,6 +282,43 @@ def make_xdmffiles(results_file, names=None):
                 logger.info(f"Could not save {name}")
 
 
+def ecg_recovery(
+    *,
+    v: dolfin.Function,
+    sigma_b: Union[dolfin.Constant, float],
+    mesh: dolfin.Mesh,
+    dx: Optional[dolfin.Measure] = None,
+    point: Optional[np.ndarray] = None,
+    r: Optional[dolfin.Function] = None,
+) -> float:
+    if dx is None:
+        dx = dolfin.dx(domain=mesh)
+    if r is None:
+        r = dolfin.SpatialCoordinate(mesh) - dolfin.Constant(point)
+    r3 = ufl.sqrt((r**2)) ** 3
+    return (1 / (4 * ufl.pi * sigma_b)) * dolfin.assemble(
+        (ufl.inner(ufl.grad(v), r) / r3) * dx,
+    )
+
+
+def ecg(
+    voltage: Iterable[dolfin.Function],
+    sigma_b: float,
+    point1: tuple[float, float, float],
+    point2: tuple[float, float, float],
+    mesh: Optional[dolfin.Mesh] = None,
+) -> list[float]:
+    values = []
+    for v in voltage:
+        if mesh is None:
+            mesh = v.function_space().mesh()
+        phi1 = ecg_recovery(v=v, sigma_b=sigma_b, point=point1, mesh=mesh)
+        phi2 = ecg_recovery(v=v, sigma_b=sigma_b, point=point2, mesh=mesh)
+        values.append(phi1 - phi2)
+
+    return values
+
+
 def plot_population(results, outdir, num_models, reset_time=True):
     plt.rcParams["svg.fonttype"] = "none"
     plt.rc("axes", labelsize=13)

src/simcardems/slabgeometry.py

@@ -38,12 +38,12 @@ class SlabGeometry(BaseGeometry):
     @staticmethod
     def default_markers() -> Dict[str, Tuple[int, int]]:
         return {
-            "X0": (2, 1),
+            "X0": (1, 2),
             "X1": (2, 2),
-            "Y0": (2, 3),
-            "Y1": (2, 4),
-            "Z0": (2, 5),
-            "Z1": (2, 6),
+            "Y0": (3, 2),
+            "Y1": (4, 2),
+            "Z0": (5, 2),
+            "Z1": (6, 2),
         }
 
     def _default_microstructure(
@@ -128,14 +128,14 @@ def create_slab_facet_function(
     ffun = dolfin.MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
     ffun.set_all(0)
 
-    x0.mark(ffun, markers["X0"][1])
-    x1.mark(ffun, markers["X1"][1])
+    x0.mark(ffun, markers["X0"][0])
+    x1.mark(ffun, markers["X1"][0])
 
-    y0.mark(ffun, markers["Y0"][1])
-    y1.mark(ffun, markers["Y1"][1])
+    y0.mark(ffun, markers["Y0"][0])
+    y1.mark(ffun, markers["Y1"][0])
 
-    z0.mark(ffun, markers["Z0"][1])
-    z1.mark(ffun, markers["Z1"][1])
+    z0.mark(ffun, markers["Z0"][0])
+    z1.mark(ffun, markers["Z1"][0])
     return ffun