Aliev-Panfilov cardiac tissue
The Aliev-Panfilov model of cardiac tissue, started from a broken wave that curls into a spiral.
- Name:
aliev_panfilov - Field:
biology/physiology - Formalism: ordinary differential equations
- References: Aliev and Panfilov 1996; Datseris and Parlitz 2022
Factory
@register_model(
field='biology/physiology',
tags=('field', 'excitable media', 'spirals', 'cardiac', 'pde'),
references=('Aliev and Panfilov 1996', 'Datseris and Parlitz 2022'),
defaults={'until': 200.0, 'dt': 2.0, 'solver': RK4Solver(dt=0.05)},
)
def aliev_panfilov(
shape: Sequence[int] = (100, 100),
*,
D: float = 1.0,
k: float = 8.0,
a: float = 0.15,
epsilon0: float = 0.002,
mu1: float = 0.2,
mu2: float = 0.3,
) -> System:
"""Return the Aliev-Panfilov model of cardiac tissue, started from a broken wave that curls into a spiral.
A spiral wave in heart muscle is a re-entrant arrhythmia, a tachycardia, and its breakup into many waves models
fibrillation. The upper half of the tissue starts excited and its left half refractory, so the free end of the
wave rotates.
Args:
shape:
The number of rows and columns of the lattice with zero-flux edges.
D:
The diffusion of the potential.
k:
The strength of the excitation.
a:
The excitation threshold.
epsilon0:
The recovery rate far from the excited state.
mu1:
The growth of the recovery rate with the recovery.
mu2:
The offset of the recovery rate.
Returns:
system:
The model, with the potential `tissue.u` and the recovery `tissue.v`.
Examples:
>>> from skcomplex.models.biology.physiology import aliev_panfilov
>>> aliev_panfilov((20, 20)).elements[0].values['k']
8.0
"""
rows, columns = shape
u, v = np.zeros(shape), np.zeros(shape)
u[: rows // 2, :] = 1.0
v[:, : columns // 2] = 2.0
values = {'u': u, 'v': v, 'k': k, 'a': a, 'epsilon0': epsilon0, 'mu1': mu1, 'mu2': mu2}
space = LatticeSpace((rows, columns), boundary='reflecting')
tissue = FieldElement('tissue', element=AlievPanfilovElement, space=space, values=values)
diffusion = DiffusionInteraction('diffusion', target='tissue', state='u', coefficient=D)
return System(elements=[tissue], interactions=[diffusion])
Components
class AlievPanfilovElement(Element):
"""Cardiac tissue of Aliev and Panfilov, an excitable medium with the restitution of heart muscle.
The potential follows `u' = -k u (u - a) (u - 1) - u v` and the recovery `v' = (epsilon0 + mu1 v / (u + mu2))
(-v - k u (u - a - 1))`, and a diffusion of the potential couples the cells of the tissue. The recovery rate
depends on the state, so the action potential lasts longer than in the FitzHugh-Nagumo medium and its duration
adapts to the pacing rate.
Args:
label:
The label of the tissue.
visual:
How the element is drawn, or `None` for `default_visual_`.
u:
The initial potential.
v:
The initial recovery.
k:
The strength of the excitation.
a:
The excitation threshold.
epsilon0:
The recovery rate far from the excited state.
mu1:
The growth of the recovery rate with the recovery.
mu2:
The offset of the recovery rate.
Examples:
>>> from skcomplex.models.biology.physiology import AlievPanfilovElement
>>> AlievPanfilovElement('tissue').k
8.0
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
u: Annotated[float, StateSpec()] = 0.0,
v: Annotated[float, StateSpec()] = 0.0,
k: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 8.0,
a: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.15,
epsilon0: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.002,
mu1: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.2,
mu2: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.3,
) -> None:
super().__init__(label, visual=visual)
self.u = u
self.v = v
self.k = k
self.a = a
self.epsilon0 = epsilon0
self.mu1 = mu1
self.mu2 = mu2
def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return the diffusing potential and its state-dependent recovery.
Args:
s:
The namespace of the symbols the laws can use.
Returns:
laws:
The laws of the component.
"""
recovery_rate = s.epsilon0 + s.mu1 * s.v / (s.u + s.mu2)
return [
RateLaw(s.u, -s.k * s.u * (s.u - s.a) * (s.u - 1) - s.u * s.v),
RateLaw(s.v, recovery_rate * (-s.v - s.k * s.u * (s.u - s.a - 1))),
]
class DiffusionInteraction(ConservationInteraction):
"""Diffusion `du/dt = D laplacian(u)` of a state `u` of a field, the conservation law of the flux `-D gradient(u)`.
The coefficient `D` is uniform in space, a number or a function of the namespace that returns an expression of
parameters and time. The interaction lowers to the Laplacian of the field, which its discretization evaluates, so
the total of the state on a periodic field stays constant and a point mass spreads as the heat kernel.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the field whose state diffuses.
state:
The name of the state of the field that diffuses.
coefficient:
A function of the namespace that returns the diffusion coefficient, or a number.
Raises:
ValidationError:
The state is not a name, or the coefficient is neither callable nor a number.
Examples:
>>> from typing import Annotated
>>> import numpy as np
>>> from skcomplex.core.elements import Element, FieldElement
>>> from skcomplex.core.interactions import DiffusionInteraction
>>> from skcomplex.core.spaces import LatticeSpace
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.specs import StateSpec
>>> class CellElement(Element):
... def __init__(self, label, u: Annotated[float, StateSpec()] = 0.0):
... super().__init__(label)
... self.u = u
>>> space = LatticeSpace((40,), spacing=0.1, boundary='periodic')
>>> field = FieldElement('rod', element=CellElement, space=space, values={'u': np.eye(40)[20]})
>>> heat = DiffusionInteraction('heat', target='rod', state='u', coefficient=0.1)
>>> trajectory = Simulation(System(elements=[field], interactions=[heat]), dt=0.01).check().prepare()
>>> round(float(trajectory.run(until=1.0).trajectory_.values('rod.u')[-1].sum()), 6)
1.0
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
state: str = 'u',
coefficient: Callable[[SimpleNamespace], object] | float = 1.0,
) -> None:
super().__init__(label, target=target, state=state)
self.coefficient = coefficient
def _check(self) -> None:
"""Raise an error for a state that is not a name, or a coefficient that is neither callable nor a number.
Raises:
ValidationError:
The type declares states, the state is not a name, or the coefficient is neither callable nor a number.
"""
super()._check()
if not (callable(self.coefficient) or isinstance(self.coefficient, int | float)):
error_msg = f'The coefficient of {self.label!r} is {self.coefficient!r}. Give a function or a number.'
raise ValidationError(error_msg)
def _get_transport_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the rate `D laplacian(u)` of the state, the divergence of the diffusive flux.
Args:
namespace:
The namespace of the symbols the laws can use, with the namespace of the field.
Returns:
laws:
The rate law of the state.
"""
coefficient = self.coefficient(namespace) if callable(self.coefficient) else self.coefficient
state = getattr(self._get_member(namespace, 'target'), self.state)
return [RateLaw(state, coefficient * laplacian(state))]
Equations
tissue.u[i]: diff(u_tissue(t), t) == (a_tissue - u_tissue(t))*k_tissue*(u_tissue(t) - 1)*u_tissue(t) - u_tissue(t)*v_tissue(t) + 1.0*laplacian(u_tissue(t))
tissue.v[i]: diff(v_tissue(t), t) == ((a_tissue - u_tissue(t) + 1)*k_tissue*u_tissue(t) - v_tissue(t))*(epsilon0_tissue + mu1_tissue*v_tissue(t)/(mu2_tissue + u_tissue(t)))