FitzHugh-Nagumo spirals
The FitzHugh-Nagumo excitable medium, started from a broken wave that curls into a spiral.
- Name:
fitzhugh_nagumo - Field:
biology/physiology - Formalism: ordinary differential equations
- References: FitzHugh 1961; Nagumo, Arimoto, and Yoshizawa 1962; Gros 2015
Factory
@register_model(
field='biology/physiology',
tags=('field', 'excitable media', 'spirals', 'pde'),
references=('FitzHugh 1961', 'Nagumo, Arimoto, and Yoshizawa 1962', 'Gros 2015'),
defaults={'until': 200.0, 'dt': 2.0, 'solver': RK4Solver(dt=0.05)},
)
def fitzhugh_nagumo(
shape: Sequence[int] = (100, 100),
*,
D: float = 1.0,
a: float = 0.7,
b: float = 0.8,
epsilon: float = 0.02,
) -> System:
"""Return the FitzHugh-Nagumo excitable medium, started from a broken wave that curls into a spiral.
The upper half starts excited and the 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 `u`.
a:
The offset of the recovery.
b:
The damping of the recovery.
epsilon:
The time scale of the recovery.
Returns:
system:
The model, with the potential `medium.u` and the recovery `medium.v`.
Examples:
>>> from skcomplex.models.biology.physiology import fitzhugh_nagumo
>>> fitzhugh_nagumo((20, 20)).interactions[0].epsilon
0.02
"""
rows, columns = shape
u, v = np.full(shape, -1.2), np.full(shape, -0.62)
u[: rows // 2, :] = 1.5
v[:, : columns // 2] = 1.0
space = LatticeSpace((rows, columns), boundary='reflecting')
medium = FieldElement('medium', element=FitzHughNagumoElement, space=space, values={'u': u, 'v': v})
kinetics = FitzHughNagumoInteraction('kinetics', target='medium', fast='u', slow='v', a=a, b=b, epsilon=epsilon)
diffusion = DiffusionInteraction('diffusion', target='medium', state='u', coefficient=D)
return System(elements=[medium], interactions=[kinetics, diffusion])
Components
class FitzHughNagumoElement(Element):
"""Cell of an excitable medium, with a potential `u` and a recovery `v`.
The cell has no laws of its own. FitzHugh-Nagumo kinetics and a diffusion of the potential change it.
Args:
label:
The label of the cell.
visual:
How the element is drawn, or `None` for `default_visual_`.
u:
The initial potential.
v:
The initial recovery.
Examples:
>>> from skcomplex.models.biology.physiology import FitzHughNagumoElement
>>> FitzHughNagumoElement('cell', u=1.5).u
1.5
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
u: Annotated[float, StateSpec()] = 0.0,
v: Annotated[float, StateSpec()] = 0.0,
) -> None:
super().__init__(label, visual=visual)
self.u = u
self.v = v
class FitzHughNagumoInteraction(Interaction):
"""FitzHugh-Nagumo kinetics of a fast variable `v` and a slow recovery `w` of a target.
The fast variable follows `v' = v - v³/3 - w + I` and the recovery `w' = epsilon (v + a - b w)`. On a field the
kinetics act on every cell, and a diffusion of the fast variable makes the medium excitable.
Args:
label:
The label of the interaction.
target:
The path of the element or field whose states follow the kinetics.
fast:
The name of the fast state, such as a membrane potential.
slow:
The name of the slow recovery state.
a:
The offset of the recovery.
b:
The damping of the recovery.
epsilon:
The time scale of the recovery.
I:
The applied current.
Raises:
ValidationError:
A state name is not a name, or the fast and the slow state are the same.
Examples:
>>> from typing import Annotated
>>> from skcomplex.core.elements import Element
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.models.biology.common.excitable import FitzHughNagumoInteraction
>>> from skcomplex.specs import StateSpec
>>> class CellElement(Element):
... def __init__(
... self, label, visual=None, *, v: Annotated[float, StateSpec()] = -1.2,
... w: Annotated[float, StateSpec()] = -0.6,
... ):
... super().__init__(label, visual=visual)
... self.v = v
... self.w = w
>>> kinetics = FitzHughNagumoInteraction('kinetics', target='cell', I=0.5)
>>> system = System(elements=[CellElement('cell')], interactions=[kinetics])
>>> trajectory = Simulation(system, dt=0.1).check().prepare().run(until=100.0).trajectory_
>>> bool(trajectory.values('cell.v').max() > 1.5)
True
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
fast: str = 'v',
slow: str = 'w',
a: Annotated[float, ParameterSpec()] = 0.7,
b: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.8,
epsilon: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.08,
I: Annotated[float, ParameterSpec()] = 0.0,
) -> None:
super().__init__(label)
self.target = target
self.fast = fast
self.slow = slow
self.a = a
self.b = b
self.epsilon = epsilon
self.I = I
def _check(self) -> None:
"""Raise an error for state names that are not names or that coincide.
Raises:
ValidationError:
A state name is not a name, or the fast and the slow state are the same.
"""
super()._check()
names = (self.fast, self.slow)
if not all(isinstance(name, str) and name.isidentifier() for name in names) or self.fast == self.slow:
error_msg = f'{self.label!r} reads the states {names!r}. Give two different names of states.'
raise ValidationError(error_msg)
def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return the cubic fast kinetics and the slow linear recovery.
Args:
s:
The namespace of the symbols the laws can use, with the namespace of the target.
Returns:
laws:
The rate laws of the fast and the slow state.
"""
v, w = getattr(s.target, self.fast), getattr(s.target, self.slow)
return [RateLaw(v, v - v**3 / 3 - w + s.I), RateLaw(w, s.epsilon * (v + s.a - s.b * w))]
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
medium.u[i]: diff(u_medium(t), t) == -1/3*u_medium(t)^3 + I_kinetics + 1.0*laplacian(u_medium(t)) + 1.0*u_medium(t) - v_medium(t)
medium.v[i]: diff(v_medium(t), t) == -(b_kinetics*v_medium(t) - a_kinetics - u_medium(t))*epsilon_kinetics