S1S2 protocol on tissue#
from typing import NamedTuple
from pathlib import Path
import beat.single_cell
import dolfin
import matplotlib.pyplot as plt
import numpy as np
import gotranx
import beat
results_folder = Path("results-s1s2-tissue")
save_every_ms = 1.0
dimension = 2
transverse = False
end_time = 20.0
dt = 0.05
overwrite = False
stim_amp = 500.0
mesh_unit = "cm"
# Load mesh
mesh_unit = mesh_unit
dx = 0.02 * beat.units.ureg("cm").to(mesh_unit).magnitude
L = 1.0 * beat.units.ureg("cm").to(mesh_unit).magnitude
data = beat.geometry.get_2D_slab_geometry(Lx=L, Ly=L, dx=dx)
V = dolfin.FunctionSpace(data.mesh, "CG", 1)
save_freq = round(save_every_ms / dt)
print("Running model")
Running model
# Load the model
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
model = tentusscher_panfilov_2006_epi_cell.__dict__
# Surface to volume ratio
chi = 1400.0 * beat.units.ureg("cm**-1")
# Membrane capacitance
C_m = 1.0 * beat.units.ureg("uF/cm**2")
fun = model["forward_generalized_rush_larsen"]
y = model["init_state_values"]()
time = dolfin.Constant(0.0)
CI1 = 350.0
CI0 = 50.0
CI = dolfin.Constant(CI1)
CIincr = 50.0
parameters = model["init_parameter_values"](stim_amplitude=0.0)
subdomain_data = dolfin.MeshFunction(
"size_t", data.mesh, data.mesh.topology().dim() - 1
)
subdomain_data.set_all(0)
marker = 1
dolfin.CompiledSubDomain("x[0] < 2 * dx", dx=dx).mark(subdomain_data, 1)
I_s = beat.stimulation.define_stimulus(
mesh=data.mesh,
chi=chi,
time=time,
subdomain_data=subdomain_data,
marker=marker,
amplitude=stim_amp,
mesh_unit=mesh_unit,
PCL=CI,
start=5.0,
duration=5.0,
)
V_ode = dolfin.FunctionSpace(data.mesh, "Lagrange", 1)
parameters_ode = np.zeros((len(parameters), V_ode.dim()))
parameters_ode.T[:] = parameters
M = beat.conductivities.define_conductivity_tensor(
f0=data.f0,
**beat.conductivities.default_conductivities("Niederer"),
)
pde = beat.MonodomainModel(time=time, mesh=data.mesh, M=M, I_s=I_s, C_m=C_m.magnitude)
ode = beat.odesolver.DolfinODESolver(
v_ode=dolfin.Function(V_ode),
v_pde=pde.state,
fun=fun,
init_states=y,
parameters=parameters_ode,
num_states=len(y),
v_index=model["state_index"]("V"),
)
solver = beat.MonodomainSplittingSolver(pde=pde, ode=ode)
Cannot determine geometric dimension from expression.
fname = (results_folder / "V.xdmf").as_posix()
beat.postprocess.delete_old_file(fname)
# Pick 12 uniformly distributed points on the mesh
points = np.array(
[
[0.1, 0.1],
[0.1, 0.5],
[0.1, 0.9],
[0.5, 0.1],
[0.5, 0.5],
[0.5, 0.9],
[0.9, 0.1],
[0.9, 0.5],
[0.9, 0.9],
]
)
V_values: list[list[float]] = [[] for _ in points]
times = []
stim_values = []
def save(t):
v = solver.pde.state.vector().get_local()
print(f"Solve for {t=:.2f}, {v.max() =}, {v.min() = }")
with dolfin.XDMFFile(data.mesh.mpi_comm(), fname) as xdmf:
xdmf.parameters["functions_share_mesh"] = True
xdmf.parameters["rewrite_function_mesh"] = False
xdmf.write_checkpoint(
solver.pde.state,
"V",
float(t),
dolfin.XDMFFile.Encoding.HDF5,
True,
)
times.append(t)
stim_values.append(I_s.expr(0))
for i, p in enumerate(points):
V_values[i].append(solver.pde.state(p))
fig, ax = plt.subplots(2, 1)
ax[0].plot(times, stim_values)
for p in points:
ax[1].plot(times, V_values[i], label=f"{p}")
ax[1].set_xlabel("Time [ms]")
ax[1].set_ylabel("V [mV]")
ax[1].legend(title="Point")
ax[0].set_title("Stimulus")
fig.savefig(results_folder / "V.png")
plt.close()
print("Saved to ", results_folder / "V.png")
np.save(
results_folder / "V.npy",
{"V": V_values, "time": times, "stim": stim_values},
allow_pickle=True,
)
save_freq = int(1.0 / dt)
num_beats = 5
CIs = np.arange(CI1, CI0 - CIincr, -CIincr)
def run():
t = 0.0
i = 0
for CI_value in CIs:
CI.assign(CI_value)
for _ in range(num_beats):
for ti in np.arange(0, CI_value, dt):
if i % save_freq == 0:
save(t)
solver.step((t, t + dt))
i += 1
t += dt
time.assign(ti)
if t > end_time:
print("Terminating")
return
run()
Solve for t=0.00, v.max() =-85.22999999999996, v.min() = -85.23000000000009
Saved to results-s1s2-tissue/V.png
Calling FFC just-in-time (JIT) compiler, this may take some time.
Solve for t=1.00, v.max() =-85.23053249836725, v.min() = -85.23053249837004
Saved to results-s1s2-tissue/V.png
Solve for t=2.00, v.max() =-85.23135471702182, v.min() = -85.23135471702643
Saved to results-s1s2-tissue/V.png
Solve for t=3.00, v.max() =-85.23241172050528, v.min() = -85.2324117205114
Saved to results-s1s2-tissue/V.png
Solve for t=4.00, v.max() =-85.23365847677802, v.min() = -85.23365847678531
Saved to results-s1s2-tissue/V.png
Solve for t=5.00, v.max() =-85.23505794261004, v.min() = -85.23505794261837
Saved to results-s1s2-tissue/V.png
Solve for t=6.00, v.max() =-52.91202598676395, v.min() = -85.33561627986877
Saved to results-s1s2-tissue/V.png
Solve for t=7.00, v.max() =18.511168737914076, v.min() = -85.24776629984419
Saved to results-s1s2-tissue/V.png
Solve for t=8.00, v.max() =49.17545943437013, v.min() = -86.66689834888845
Saved to results-s1s2-tissue/V.png
Solve for t=9.00, v.max() =50.63883486161976, v.min() = -87.86953943548767
Saved to results-s1s2-tissue/V.png
Solve for t=10.00, v.max() =52.34756889091984, v.min() = -88.5354778364135
Saved to results-s1s2-tissue/V.png
Solve for t=11.00, v.max() =27.86854032310785, v.min() = -88.65149435183953
Saved to results-s1s2-tissue/V.png
Solve for t=12.00, v.max() =27.23735538766364, v.min() = -88.697414639174
Saved to results-s1s2-tissue/V.png
Solve for t=13.00, v.max() =27.25342949764352, v.min() = -88.63419317764459
Saved to results-s1s2-tissue/V.png
Solve for t=14.00, v.max() =27.223252441060097, v.min() = -88.83500927786102
Saved to results-s1s2-tissue/V.png
Solve for t=15.00, v.max() =27.377706512871036, v.min() = -88.78543887018667
Saved to results-s1s2-tissue/V.png
Solve for t=16.00, v.max() =27.347179025017653, v.min() = -88.63893378695285
Saved to results-s1s2-tissue/V.png
Solve for t=17.00, v.max() =27.208558381352358, v.min() = -88.80058156972596
Saved to results-s1s2-tissue/V.png
Solve for t=18.00, v.max() =27.481368467805247, v.min() = -88.81855282644287
Saved to results-s1s2-tissue/V.png
Solve for t=19.00, v.max() =27.39898270324071, v.min() = -88.67719991772334
Saved to results-s1s2-tissue/V.png
Terminating