Barkley excitable medium with spirals and breakup
Barkley's excitable medium, started from orthogonal stripes whose crossing curls into a spiral wave.
- Name:
barkley - Field:
biology/physiology - Formalism: ordinary differential equations
- References: Barkley 1991; Datseris and Parlitz 2022
Select this system
from skcomplex.systems import create_simulation, create_system, describe
system = create_system('barkley')
entry = describe('barkley')
trajectory = create_simulation('barkley').prepare().run().trajectory_
Factory
@register_system(
field='biology/physiology',
tags=('field', 'excitable media', 'spirals', 'spatiotemporal chaos', 'pde'),
references=('Barkley 1991', 'Datseris and Parlitz 2022'),
defaults={'until': 400.0, 'dt': 4.0, 'solver': RK4Solver(dt=0.2)},
)
def barkley(
shape: Sequence[int] = (64, 64),
*,
spacing: float = 3.0,
d: float = 1.0,
a: float = 0.6,
b: float = 0.05,
epsilon: float = 0.02,
gamma: float = 1.0,
) -> System:
"""Return Barkley's excitable medium, started from orthogonal stripes whose crossing curls into a spiral wave.
The excitation `u` and the recovery `w` start as orthogonal stripes of zeros and ones, so the excited stripe runs
into recovering tissue on one side only and its free end rotates. With `(a, b, d, epsilon) = (0.6, 0.05, 1, 0.02)`
the spiral is stable for `gamma = 1`. With `(0.7, 0.05, 1, 0.08)` and `gamma = 3` it breaks up into spatiotemporal
chaos. The source runs a domain of 300 by 300, and the default lattice covers 192 by 192 with a spacing of 3, wide
enough for the stable spiral.
Args:
shape:
The number of rows and columns of the lattice with zero-flux edges.
spacing:
The distance between neighbouring sites.
d:
The diffusion of the excitation.
a:
The inverse slope of the threshold.
b:
The offset of the threshold.
epsilon:
The ratio of the time scales of the recovery and the excitation.
gamma:
The exponent of the excitation in the recovery.
Returns:
system:
The system, with the excitation `medium.u` and the recovery `medium.w`.
Examples:
>>> from skcomplex.systems.biology.physiology import barkley
>>> barkley((20, 20), gamma=3.0).elements[0].values['gamma']
3.0
"""
rows, columns = shape
u, w = np.zeros(shape), np.zeros(shape)
u[:, : columns // 2] = 1.0
w[: rows // 2, :] = 1.0
values = {'u': u, 'w': w, 'a': a, 'b': b, 'epsilon': epsilon, 'gamma': gamma}
space = LatticeSpace((rows, columns), boundary='reflecting', spacing=spacing)
medium = FieldElement('medium', element=BarkleyElement, space=space, values=values)
diffusion = DiffusionInteraction('diffusion', target='medium', state='u', coefficient=d)
return System(elements=[medium], interactions=[diffusion])
Components
Elements:
BarkleyElement, fromElement
Interactions:
DiffusionInteraction
class BarkleyElement(Element):
"""Excitable medium of Barkley, a fast excitation `u` and a slow recovery `w` with a polynomial nullcline.
The excitation follows `u' = u (1 - u) (u - (w + b) / a)` and the recovery `w' = epsilon (u^gamma - w)`. A diffusion
of the excitation couples the sites of the medium. The threshold `(w + b) / a` rises with the recovery, so an
excited site recovers before it can fire again.
Args:
label:
The label of the medium.
visual:
How the element is drawn, or `None` for `default_visual_`.
u:
The initial excitation.
w:
The initial recovery.
a:
The inverse slope of the threshold.
b:
The offset of the threshold.
epsilon:
The ratio of the time scales of the recovery and the excitation.
gamma:
The exponent of the excitation in the recovery, 1 for Barkley's model and 3 for the variant of Bär and
Eiswirth with spiral breakup.
Examples:
>>> from skcomplex.systems.biology.physiology import BarkleyElement
>>> BarkleyElement('medium').a
0.6
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
u: Annotated[float, StateSpec()] = 0.0,
w: Annotated[float, StateSpec()] = 0.0,
a: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.6,
b: Annotated[float, ParameterSpec()] = 0.05,
epsilon: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.02,
gamma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
super().__init__(label, visual=visual)
self.u = u
self.w = w
self.a = a
self.b = b
self.epsilon = epsilon
self.gamma = gamma
def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return the cubic excitation and its slow recovery.
Args:
s:
The namespace of the symbols the laws can use.
Returns:
laws:
The laws of the component.
"""
return [
RateLaw(s.u, s.u * (1 - s.u) * (s.u - (s.w + s.b) / s.a)),
RateLaw(s.w, s.epsilon * (s.u**s.gamma - s.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.components.elements import Element, FieldElement
>>> from skcomplex.core.components.interactions import DiffusionInteraction
>>> from skcomplex.core.components.spaces import LatticeSpace
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.runtime import Simulation
>>> from skcomplex.core.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).prepare()
>>> round(float(trajectory.run(until=1.0).trajectory_.values('rod.u')[-1].sum()), 6)
1.0
"""
law_classes: ClassVar[tuple[type[BaseLaw], ...]] = (FluxLaw, RateLaw)
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec(elements=(ElementKind.ANY,), states=('state',))],
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) == ((b_medium + w_medium(t))/a_medium - u_medium(t))*(u_medium(t) - 1)*u_medium(t) + 1.0*laplacian(u_medium(t))
medium.w[i]: diff(w_medium(t), t) == epsilon_medium*(u_medium(t)^gamma_medium - w_medium(t))