Niederer benchmark

Niederer benchmark#

In this example we will use the same setup as in the Niederer benchmark [LGA+15].

from pathlib import Path
import json

import beat.conductivities
import dolfin
import numpy.typing as npt
import pyvista
import gotranx
import matplotlib.pyplot as plt

import beat
import beat.viz

def setup_initial_conditions() -> npt.NDArray:
    ic = {
        "V": -85.23,  # mV
        "Xr1": 0.00621,
        "Xr2": 0.4712,
        "Xs": 0.0095,
        "m": 0.00172,
        "h": 0.7444,
        "j": 0.7045,
        "d": 3.373e-05,
        "f": 0.7888,
        "f2": 0.9755,
        "fCass": 0.9953,
        "s": 0.999998,
        "r": 2.42e-08,
        "Ca_i": 0.000126,  # millimolar
        "R_prime": 0.9073,
        "Ca_SR": 3.64,  # millimolar
        "Ca_ss": 0.00036,  # millimolar
        "Na_i": 8.604,  # millimolar
        "K_i": 136.89,  # millimolar
    }
    values = model.init_state_values(**ic)
    return values

dx = 0.5
dt = 0.05
# Increase T to 100 to reproduce Niederer benchmark
T = 20.0

here = Path.cwd()

model_path = Path("tentusscher_panfilov_2006_epi_cell.py")
if not model_path.is_file():
    here = Path.cwd()
    ode = gotranx.load_ode(
        here
        / ".."
        / "odes"
        / "tentusscher_panfilov_2006"
        / "tentusscher_panfilov_2006_epi_cell.ode"
    )
    code = gotranx.cli.gotran2py.get_code(
        ode, scheme=[gotranx.schemes.Scheme.forward_generalized_rush_larsen]
    )
    model_path.write_text(code)

import tentusscher_panfilov_2006_epi_cell as model

fun = model.forward_generalized_rush_larsen
init_states = setup_initial_conditions()
parameters = model.init_parameter_values(stim_amplitude=0.0)

mesh_unit = "mm"

Lx = 20.0 * beat.units.ureg("mm").to(mesh_unit).magnitude
Ly = 7 * beat.units.ureg("mm").to(mesh_unit).magnitude
Lz = 3 * beat.units.ureg("mm").to(mesh_unit).magnitude
dx_mm = dx * beat.units.ureg("mm").to(mesh_unit).magnitude

geo = beat.geometry.get_3D_slab_geometry(
    Lx=Lx,
    Ly=Ly,
    Lz=Lz,
    dx=dx_mm,
)
ode_space = dolfin.FunctionSpace(geo.mesh, "Lagrange", 1)

pyvista.start_xvfb()
plotter = pyvista.Plotter()
topology, cell_types, x = beat.viz.create_vtk_structures(ode_space)
grid = pyvista.UnstructuredGrid(topology, cell_types, x)
plotter.add_mesh(grid, show_edges=True)
plotter.show_grid()
plotter.add_axes(line_width=5)
plotter.show_axes()
plotter.view_xy()

if not pyvista.OFF_SCREEN:
    plotter.show()
else:
    figure = plotter.screenshot("niederer_mesh.png")

# Surface to volume ratio
conductivities = beat.conductivities.default_conductivities("Niederer")
# # Membrane capacitance
C_m = 1.0 * beat.units.ureg("uF/cm**2")

time = dolfin.Constant(0.0)
L = 1.5 * beat.units.ureg("mm").to(mesh_unit).magnitude
S1_marker = 1
S1_subdomain = dolfin.CompiledSubDomain(
    "x[0] <= L + DOLFIN_EPS && x[1] <= L + DOLFIN_EPS && x[2] <= L + DOLFIN_EPS",
    L=L,
)
S1_markers = dolfin.MeshFunction("size_t", geo.mesh, geo.mesh.topology().dim())
S1_subdomain.mark(S1_markers, S1_marker)

I_s = beat.stimulation.define_stimulus(
    mesh=geo.mesh,
    chi=conductivities["chi"],
    time=time,
    subdomain_data=S1_markers,
    marker=S1_marker,
    mesh_unit=mesh_unit,
    amplitude=50_000.0,
)


M = beat.conductivities.define_conductivity_tensor(
    f0=geo.f0,
    **conductivities,
)

params = {"preconditioner": "sor", "use_custom_preconditioner": False}

pde = beat.MonodomainModel(
    time=time,
    mesh=geo.mesh,
    M=M,
    I_s=I_s,
    params=params,
    C_m=C_m.to(f"uF/{mesh_unit}**2").magnitude,
)
ode = beat.odesolver.DolfinODESolver(
    v_ode=dolfin.Function(ode_space),
    v_pde=pde.state,
    fun=fun,
    init_states=init_states,
    parameters=parameters,
    num_states=len(init_states),
    v_index=model.state_index("V"),
)

Cannot determine geometric dimension from expression.

Calling FFC just-in-time (JIT) compiler, this may take some time.

solver = beat.MonodomainSplittingSolver(pde=pde, ode=ode, theta=0.0)
output_dir = Path("results-niederer-benchmark")
output_dir.mkdir(exist_ok=True)
filename = output_dir / f"results-{dt}-{dx}.xdmf"
filename.unlink(missing_ok=True)
filename.with_suffix(".h5").unlink(missing_ok=True)

xdmf = dolfin.XDMFFile(geo.mesh.mpi_comm(), filename.as_posix())
xdmf.parameters["functions_share_mesh"] = True
xdmf.parameters["rewrite_function_mesh"] = False
xdmf.write(geo.mesh)

points = {
    "P1": (0, 0, 0),
    "P2": (0.0, Ly, 0.0),
    "P3": (Lx, 0.0, 0.0),
    "P4": (Lx, Ly, 0.0),
    "P5": (0.0, 0.0, Lz),
    "P6": (0.0, Ly, Lz),
    "P7": (Lx, 0.0, Lz),
    "P8": (Lx, Ly, Lz),
    "P9": (Lx / 2, Ly / 2, Lz / 2),
}

activation_times = {p: -1.0 for p in points}
save_freq = int(1.0 / dt)
i = 0
plotter_voltage = pyvista.Plotter()
viridis = plt.get_cmap("viridis")
grid.point_data["V"] = solver.pde.state.vector().get_local()
grid.set_active_scalars("V")
renderer = plotter_voltage.add_mesh(
    grid,
    show_edges=True,
    lighting=False,
    cmap=viridis,
    clim=[-90.0, 40.0],
)
gif_file = Path("niederer_benchmark.gif")
gif_file.unlink(missing_ok=True)
plotter_voltage.open_gif(gif_file.as_posix())

T = 20
# T = 100  # Change to 100 to reproduce Niederer benchmark
t = 0.0
while t < T + 1e-12 and any(at < 0.0 for at in activation_times.values()):
    v = solver.pde.state.vector().get_local()
    if i % save_freq == 0:
        print(f"Solve for {t=:.2f}, {v.max() =}, {v.min() = }")
        print(activation_times)
        xdmf.write_checkpoint(
            solver.pde.state,
            "V",
            float(t),
            dolfin.XDMFFile.Encoding.HDF5,
            True,
        )
        grid.point_data["V"] = solver.pde.state.vector().get_local()
        plotter_voltage.write_frame()
    solver.step((t, t + dt))

    for p in points:
        value = beat.utils.peval(solver.pde.state, points[p])
        if value > 0.0 and activation_times[p] < 0.0:
            activation_times[p] = t
    i += 1
    t += dt

plotter_voltage.close()

Solve for t=0.00, v.max() =-85.22999999999993, v.min() = -85.23000000000002
{'P1': -1.0, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}
Calling FFC just-in-time (JIT) compiler, this may take some time.

Solve for t=1.00, v.max() =-16.870378443055664, v.min() = -88.13697024957764
{'P1': -1.0, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=2.00, v.max() =60.777926916946484, v.min() = -88.81234270816279
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=3.00, v.max() =40.468888963769324, v.min() = -89.58747927438704
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=4.00, v.max() =29.746407967798497, v.min() = -89.6450675674948
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=5.00, v.max() =32.17405983997671, v.min() = -89.46845566962915
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=6.00, v.max() =30.909440035882014, v.min() = -88.75412142316867
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=7.00, v.max() =30.319493378426568, v.min() = -89.03863299063237
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=8.00, v.max() =34.18807603301506, v.min() = -88.92247667576746
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=9.00, v.max() =29.10333807548235, v.min() = -89.43365419802898
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=10.00, v.max() =31.951746638750038, v.min() = -89.48251977542454
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=11.00, v.max() =30.343642002885034, v.min() = -89.5574792354038
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=12.00, v.max() =27.72546844255657, v.min() = -89.9414863143836
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=13.00, v.max() =37.28698160430252, v.min() = -90.00985538164105
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': -1.0, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=14.00, v.max() =38.75304874745052, v.min() = -90.68118525203792
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': 13.950000000000063, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=15.00, v.max() =36.24740836556221, v.min() = -90.8659044587285
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': 13.950000000000063, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=16.00, v.max() =40.57379334341332, v.min() = -91.31466501661228
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': 13.950000000000063, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=17.00, v.max() =36.46403848989674, v.min() = -91.31786039094672
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': 13.950000000000063, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=18.00, v.max() =35.99765749677111, v.min() = -91.60532761697284
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': 13.950000000000063, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=19.00, v.max() =36.94440206132663, v.min() = -91.96602183941026
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': 13.950000000000063, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

Solve for t=20.00, v.max() =33.35964911043149, v.min() = -92.03969706814676
{'P1': 1.2500000000000004, 'P2': -1.0, 'P3': -1.0, 'P4': -1.0, 'P5': 13.950000000000063, 'P6': -1.0, 'P7': -1.0, 'P8': -1.0, 'P9': -1.0}

_

# Save activation times
activation_times["dx"] = dx
activation_times["dt"] = dt
at_file_name = output_dir / "activation_times.json"
if at_file_name.is_file():
    all_at = json.loads(at_file_name.read_text())
else:
    all_at = []
all_at.append(activation_times)
at_file_name.write_text(json.dumps(all_at, indent=2))

213

The activation times are saved in the file output-niederer-benchmark/activation_times.json. The file contains a list of dictionaries, each dictionary contains the activation times for a specific dx and dt. Here are the activation times for the different dx and dt:

dx

dt

P1

P2

P3

P4

P5

P6

P7

P8

P9

0

0.5

0.05

1.25

51.1

34.9

58.9

14.1

49.5

34

56.65

26.05

1

0.5

0.01

1.22

50.85

33.96

58.05

13.98

49.36

33.07

55.91

25.64

2

0.5

0.005

1.215

50.775

33.825

57.96

13.97

49.345

32.945

55.825

25.595

3

0.2

0.05

1.25

29.7

32.9

40.2

9.55

30

32.95

39.9

18.9

4

0.2

0.01

1.24

29.09

31.25

38.66

9.34

29.4

31.29

38.42

18.14

5

0.2

0.005

1.235

29.015

31.05

38.475

9.315

29.32

31.08

38.235

18.045

6

0.1

0.05

1.25

26.85

33.3

40.35

8.4

27.5

33.85

40.55

18.95

7

0.1

0.01

1.23

25.64

31.46

38.08

8.03

26.24

31.94

38.21

17.95

8

0.1

0.005

1.225

25.5

31.26

37.81

7.99

26.09

31.72

37.93

17.835
[CRN98]

Marc Courtemanche, Rafael J Ramirez, and Stanley Nattel. Ionic mechanisms underlying human atrial action potential properties: insights from a mathematical model. American Journal of Physiology-Heart and Circulatory Physiology, 275(1):H301–H321, 1998.

[LGA+15]

Sander Land, Viatcheslav Gurev, Sander Arens, Christoph M Augustin, Lukas Baron, Robert Blake, Chris Bradley, Sebastian Castro, Andrew Crozier, Marco Favino, and others. Verification of cardiac mechanics software: benchmark problems and solutions for testing active and passive material behaviour. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 471(2184):20150641, 2015.

[ZLH+21]

Zhaoyang Zhang, Michael B Liu, Xiaodong Huang, Zhen Song, and Zhilin Qu. Mechanisms of premature ventricular complexes caused by qt prolongation. Biophysical journal, 120(2):352–369, 2021.