Kinouchi-Copelli excitable network
The Kinouchi-Copelli network of excitable elements, whose dynamic range peaks at the critical point.
- Name:
kinouchi_copelli - Field:
biology/neuroscience - Formalism: discrete map
- References: Kinouchi and Copelli 2006; Larremore et al. 2011
Factory
@register_model(
field='biology/neuroscience',
tags=('excitable', 'network', 'criticality', 'dynamic range', 'stochastic'),
references=('Kinouchi and Copelli 2006', 'Larremore et al. 2011'),
defaults={'until': 1000, 'dt': 1},
observables=lambda system: [
ObservableSpec('activity', group='neurons', variable='phase', reduce=_excited_fraction),
],
)
def kinouchi_copelli(
graph: nx.Graph | None = None,
*,
sigma: float = 1.0,
stimulus: float = 0.0,
n_states: int = 5,
excited: float | None = None,
random_state: int | np.random.Generator | None = 0,
) -> System:
"""Return the Kinouchi-Copelli network of excitable elements, whose dynamic range peaks at the critical point.
Every element excites each quiescent neighbour with the transmission probability `p = sigma / K`, where `K` is
the mean degree of the graph, so `sigma` is the branching ratio, the mean number of elements one excitation
excites. Below `sigma = 1` activity dies out without a stimulus, and above it activity sustains itself. The
response, the mean firing rate against the stimulus rate, spans the most decades of stimulus at `sigma = 1`.
Args:
graph:
The NetworkX graph, or `None` for an Erdős-Rényi graph of 1000 nodes with mean degree 10.
sigma:
The branching ratio.
stimulus:
The rate `r` of the external stimulus per step, which excites a quiescent element with probability
`1 - exp(-r)`.
n_states:
The number of states of the cycle, one quiescent, one excited, and `n_states - 2` refractory.
excited:
The initial fraction of excited elements, with the others quiescent, or `None` for states drawn uniformly
over the cycle.
random_state:
The seed or generator of the default graph and of the initial states.
Returns:
system:
The model, with the states `neurons.phase`, one value per node. Its standard observable `activity` is the
fraction of excited elements, the firing rate per step.
Examples:
>>> import networkx as nx
>>> from skcomplex.models.biology.neuroscience import kinouchi_copelli
>>> kinouchi_copelli(nx.random_regular_graph(4, 20, seed=0), sigma=1.0).elements[0].values['p']
0.25
"""
rng = np.random.default_rng(random_state)
graph = nx.gnp_random_graph(1000, 10 / 999, seed=int(rng.integers(2**31))) if graph is None else graph
size = graph.number_of_nodes()
degree = 2 * graph.number_of_edges() / size
phase = (
rng.integers(0, n_states, size=size).astype(float)
if excited is None
else (rng.random(size) < excited).astype(float)
)
neurons = LatticeElement(
'neurons',
element=ExcitableElement,
space=NetworkSpace(graph),
values={'phase': phase, 'p': sigma / degree, 'stimulus': stimulus},
)
excitation = TransitionInteraction(
'excitation',
target='neurons',
state='phase',
source=0,
destination=EXCITED,
probability=_get_excitation,
)
return System(elements=[neurons], interactions=[excitation, *_get_cycle(n_states)], time=DiscreteTime())
Components
class ExcitableElement(Element):
"""Cyclic excitable element of the Kinouchi-Copelli network: quiescent at 0, excited at 1, refractory up to `m - 1`.
Transitions move its phase every step. A quiescent element is excited by the external stimulus with probability
`1 - exp(-r)`, where `r` is the stimulus rate per step, or by each excited neighbour with the transmission
probability `p`, so it stays quiescent with probability `exp(-r) (1 - p)^n` when `n` neighbours are excited. An
excited or refractory element moves to the next state, and the last refractory state returns to quiescence.
Args:
label:
The label of the element or the site group.
visual:
How the element is drawn, or `None` for `default_visual_`.
phase:
The initial state, 0 for quiescent, 1 for excited, and 2 to `n_states - 1` for refractory.
p:
The probability that an excited neighbour excites the element in one step.
stimulus:
The rate `r` of the external stimulus per step.
Examples:
>>> from skcomplex.models.biology.neuroscience import ExcitableElement
>>> ExcitableElement('neuron').p
0.1
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
phase: Annotated[float, StateSpec(domain=NN)] = 0,
p: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
stimulus: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
) -> None:
super().__init__(label, visual=visual)
self.phase = phase
self.p = p
self.stimulus = stimulus
class TransitionInteraction(Interaction):
"""Move of a discrete state of the target from one value to another, or between all its values by a matrix.
A member whose state holds `source` moves to `destination` at a `rate`, as a jump in continuous time, with a
`probability` in every step of discrete time, or in the step in which a `condition` holds. With `over='neighbors'`
the rate and the probability are those of one neighbour, written with `namespace.target.other` for the neighbour:
the rates of the neighbours add up, and their probabilities combine as `1 - prod(1 - p)`, the probability that at
least one of them moves the member. With a matrix instead, the state takes the values `0` to `n - 1` and steps by
the rows of a transition matrix in discrete time, the row of its current value. An `emission` matrix then draws the
`observed` state of the target from the row that the next value of the state selects, as a hidden Markov source.
With a residence time `after`, a member moves by a probability or a condition only once it has held the source
value for `after` steps, at the earliest at the end of the step that completes them.
In the `'symbolic'` form, every function reads `namespace.target`, the states and parameters of the target, and the
time `namespace.t`. In the `'agent'` form, the target is an agent group and every function takes the agent and
the model view, `f(agent, model)`, and returns a number, a truth value, or the rows of a matrix, and a rate `r`
moves an agent with probability `1 - exp(-r)` in a step. With `over='cell'`, every member acts on the other members
of its cell: it moves each one that holds the source value with the probability, the rate, or the condition
`f(agent, other, model)` of that member `agent` and the acting member `other`, each by its own draw.
The `schedule` of the agent form orders the members. Under `'synchronous'` every member decides from the states at
the start of the step, under `'sequential'` the members act one after the other in their order, and under
`'random-sequential'` in a new random order every step, each reading the values the members before it set. Under
`'member'` every member acts in its own activation, together with the rules of its element and the other member
interactions, in the order the system lists them, so one rule of a member runs its movement, transfers,
transitions, and life cycle.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose state moves.
state:
The name of the discrete state of the target.
source:
The value the state moves from, a collection of such values, or a function of the value of the state that
returns whether the state moves from it, such as `lambda value: value >= 0`, a condition in the symbolic
form. A state that moves from several values moves to the destination from each.
destination:
The value the state moves to.
rate:
The function of the rate of the move, or `None`.
probability:
The function of the probability of the move in a step, or `None`.
condition:
The function of the condition under which the state moves in a step, or `None`.
matrix:
The function of the transition matrix, one row and one column per value of the state, or `None`.
emission:
With a matrix, the function of the emission matrix, one row per value of the state and one column per value
of the observed state, or `None`.
observed:
The name of the state of the target that the emission matrix draws, or `None`.
over:
`'member'` for a rate or a probability of the member, `'neighbors'` for that of one neighbour, or `'cell'`
for the agent form whose members act on the members of their cell.
neighborhood:
The neighbourhood of `over='neighbors'`, or `None` for the neighbours of the network or the lattice.
after:
The number of steps a member holds the source value before it can move, or `None`.
schedule:
The order of the members of the agent form, `'synchronous'`, `'sequential'`, `'random-sequential'`, or
`'member'`.
form:
`'symbolic'` for laws of the states of the target, or `'agent'` for a rule over an agent group.
phase:
The phase of the step the agent form runs in, or `None` for the default phase.
Raises:
ValidationError:
The transition has not exactly one kind, the source and the destination are equal, an emission has no
matrix or no observed state, `over` is unknown or goes with a form or a kind it cannot take, the
residence time is negative or goes with a rate, a matrix, or members acting on their cell, the form or the
schedule is unknown, the symbolic form has a schedule, or a function is not callable.
Examples:
>>> from skcomplex.core.interactions import TransitionInteraction
>>> recovery = TransitionInteraction(
... 'recovery', target='hosts', state='state', source=1, destination=0, rate=lambda s: 0.5
... )
>>> recovery.check().over
'member'
>>> agents = TransitionInteraction(
... 'recovery', target='hosts', state='state', source=1, destination=2, probability=lambda a, m: 0.1,
... after=3, form='agent'
... )
>>> type(agents.check().prepare().run(SimpleNamespace()).laws_[0]).__name__
'RuleLaw'
>>> separation = TransitionInteraction(
... 'separation', target='workers', state='employer', source=lambda value: value >= 0, destination=-1,
... probability=lambda a, m: 0.05, form='agent'
... )
>>> separation.check().destination
-1
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
state: str = 'state',
source: SourceValues = 0.0,
destination: float = 1.0,
rate: Callable[..., SupportsFloat] | None = None,
probability: Callable[..., SupportsFloat] | None = None,
condition: Callable[..., object] | None = None,
matrix: Callable[..., Sequence[Sequence[object]]] | None = None,
emission: Callable[..., Sequence[Sequence[object]]] | None = None,
observed: str | None = None,
over: str = 'member',
neighborhood: BaseNeighborhood | None = None,
after: float | None = None,
schedule: str = 'synchronous',
form: str = 'symbolic',
phase: str | None = None,
) -> None:
super().__init__(label)
self.target = target
self.state = state
self.source = source
self.destination = destination
self.rate = rate
self.probability = probability
self.condition = condition
self.matrix = matrix
self.emission = emission
self.observed = observed
self.over = over
self.neighborhood = neighborhood
self.after = after
self.schedule = schedule
self.form = form
self.phase = phase
def _check_kind(self) -> str:
"""Raise an error for a transition without exactly one callable kind, or with an emission it cannot take."""
kinds = [kind for kind in KINDS if getattr(self, kind) is not None]
if len(kinds) != 1 or not callable(getattr(self, kinds[0])):
error_msg = (
f'The transition {self.label!r} has the forms {kinds}. Give exactly one of {", ".join(KINDS)}, as a '
'function of the namespace.'
)
raise ValidationError(error_msg)
unpaired = (self.emission is None) != (self.observed is None)
if unpaired or (self.emission is not None and self.matrix is None):
error_msg = (
f'The transition {self.label!r} has the emission {self.emission!r} of the state {self.observed!r}. '
'Give both with a matrix, or neither.'
)
raise ValidationError(error_msg)
return kinds[0]
def _check_form(self) -> None:
"""Raise an error for a transition without exactly one callable kind, or with arguments its kind ignores."""
kind = self._check_kind()
if self.over not in OVERS:
error_msg = f'The transition {self.label!r} runs over {self.over!r}. Use one of {", ".join(OVERS)}.'
raise ValidationError(error_msg)
if self.form not in FORMS:
error_msg = f'The form of transition {self.label!r} is {self.form!r}. Use one of {", ".join(FORMS)}.'
raise ValidationError(error_msg)
if self.schedule not in SCHEDULES or (self.form == 'symbolic' and self.schedule != 'synchronous'):
error_msg = (
f'The schedule of transition {self.label!r} is {self.schedule!r}. Use one of {", ".join(SCHEDULES)}, '
'and the agent form for another than `synchronous`.'
)
raise ValidationError(error_msg)
if self.over == 'cell' and (self.form != 'agent' or kind == 'matrix' or self.after is not None):
error_msg = (
f'The members of transition {self.label!r} act on their cell, which only the agent form does, by a '
'rate, a probability, or a condition and without a residence time.'
)
raise ValidationError(error_msg)
if self.over == 'neighbors' and (kind not in {'rate', 'probability'} or self.form == 'agent'):
error_msg = (
f'The transition {self.label!r} runs over neighbours, which only a rate or a probability of the '
'symbolic form does. In the agent form, read the neighbours of the agent in the function.'
)
raise ValidationError(error_msg)
single = not (callable(self.source) or isinstance(self.source, Collection))
if self.matrix is None and single and self.source == self.destination:
error_msg = f'The transition {self.label!r} moves from {self.source} to itself. Give distinct values.'
raise ValidationError(error_msg)
if self.after is not None and (kind not in {'probability', 'condition'} or not self.after >= 0):
error_msg = (
f'The transition {self.label!r} has the residence time {self.after!r}. Give a number of steps at least '
'0 with a probability or a condition.'
)
raise ValidationError(error_msg)
def _check_transition(self) -> None:
"""Raise an error for arguments the transition cannot take."""
self._check_form()
def _check(self) -> None:
"""Raise an error for a transition without exactly one kind, or with arguments its kind ignores.
Raises:
ValidationError:
The transition has not exactly one kind, the source and the destination are equal, an emission has
no matrix or no observed state, `over` is unknown or goes with a condition, a matrix, or the agent
form, the residence time is negative or goes with a rate or a matrix, the form is unknown, or a
function is not callable.
"""
super()._check()
self._check_transition()
def _get_moving(self, value: Expression, source: SourceValues) -> Expression:
"""Return 1 for a member whose state holds a source value, and 0 otherwise."""
if callable(source):
return where(SR(source(value)), 1, 0)
if isinstance(source, Collection):
return sum((where(value == item, 1, 0) for item in set(source)), SR(0))
return where(value == source, 1, 0)
def _get_jump(self, value: Expression, source: SourceValues, destination: float, rate: Expression) -> JumpLaw:
"""Return the jump of the state from a source value to the destination at a rate."""
if callable(source) or isinstance(source, Collection):
return JumpLaw(rate * self._get_moving(value, source), {value: destination - value})
return JumpLaw(rate * self._get_moving(value, source), {value: destination - source})
def _get_step(self, value: Expression, source: SourceValues, destination: float, moves: Expression) -> UpdateLaw:
"""Return the update of the state that moves a member from a source value where a condition holds."""
if callable(source) or isinstance(source, Collection):
return UpdateLaw(value, where(self._get_moving(value, source) > 0, where(moves, destination, value), value))
return UpdateLaw(value, where(value == source, where(moves, destination, source), value), source=source)
def _holds_source(self, value: float) -> bool:
"""Return whether the value of the state of an agent is a source value."""
if callable(self.source):
return bool(self.source(value))
if isinstance(self.source, Collection):
return value in set(self.source)
return value == self.source
def _combine_rates(self, rate: Expression) -> Expression:
"""Return the rate of the member, or the sum of the rates of its neighbours."""
return neighbors(rate, over=self.neighborhood) if self.over == 'neighbors' else rate
def _combine_probabilities(self, probability: Expression) -> Expression:
"""Return the probability of the member, or the probability that at least one neighbour moves it."""
if self.over == 'member':
return probability
return 1 - exp(neighbors(log(1 - probability), over=self.neighborhood))
def _get_draw(self) -> Expression:
"""Return the uniform draw of the transition in a step."""
return uniform(re.sub(r'\W', '_', f'transition_{self.label}'))
def _get_chain(self, namespace: SimpleNamespace, value: Expression) -> list[BasePart]:
"""Return the transition law of the state by its matrix, with the emission law of the observed state."""
laws: list[BasePart] = []
if self.matrix is not None:
laws.append(TransitionLaw(value, self.matrix(namespace), label=self.label))
if self.emission is not None and self.observed is not None:
observed = getattr(self._get_member(namespace, 'target'), self.observed)
laws.append(EmissionLaw(observed, value, self.emission(namespace), label=self.label))
return laws
def _get_held(self, agent: AgentViewProtocol, model: ModelViewProtocol) -> int:
"""Return the number of whole steps an agent has held its value, counted in the storage of the model."""
value = float(getattr(agent, self.state))
record = model.storage.get(f'{self.label}.{agent.id}')
held = record[1] + 1 if isinstance(record, tuple) and record[0] == value else 0
model.storage[f'{self.label}.{agent.id}'] = (value, held)
return held
def _draw_row(self, row: Sequence[object], model: ModelViewProtocol) -> float:
"""Return a value drawn by the probabilities of a row of a matrix."""
weights = np.asarray(row, dtype=float)
return float(model.rng.choice(weights.size, p=weights / weights.sum()))
def _evaluate[T](
self,
function: Callable[..., T],
agent: AgentViewProtocol,
other: AgentViewProtocol | None,
model: ModelViewProtocol,
) -> T:
"""Return a function of an agent, or of an agent and the member of its cell that acts on it."""
return function(agent, model) if other is None else function(agent, other, model)
def _moves_agent(
self,
agent: AgentViewProtocol,
model: ModelViewProtocol,
other: AgentViewProtocol | None = None,
) -> bool:
"""Return whether an agent that holds the source value moves in this step, moved by another one or alone."""
if self.condition is not None:
return bool(self._evaluate(self.condition, agent, other, model))
if self.rate is not None:
probability = 1.0 - math.exp(-float(self._evaluate(self.rate, agent, other, model)))
elif self.probability is not None:
probability = float(self._evaluate(self.probability, agent, other, model))
else:
probability = 0.0
return bool(model.rng.random() < probability)
def _get_next(self, agent: AgentViewProtocol, model: ModelViewProtocol) -> float | None:
"""Return the next value of the state of an agent, or `None` when it keeps its value."""
value = float(getattr(agent, self.state))
if self.matrix is not None:
return self._draw_row(self.matrix(agent, model)[int(value)], model)
held = self._get_held(agent, model) if self.after is not None else 0
if not self._holds_source(value) or (self.after is not None and held + 1 < self.after):
return None
return float(self.destination) if self._moves_agent(agent, model) else None
def _get_moves(self, agent: AgentViewProtocol, model: ModelViewProtocol) -> list[tuple[AgentViewProtocol, float]]:
"""Return the members an agent moves in this step, itself or the members of its cell, with their values."""
if self.over != 'cell':
value = self._get_next(agent, model)
return [] if value is None else [(agent, value)]
return [
(other, float(self.destination))
for other in agent.agents_here(agent.group)
if other is not agent
and self._holds_source(float(getattr(other, self.state)))
and self._moves_agent(other, model, agent)
]
def _apply(self, moves: list[tuple[AgentViewProtocol, float]], model: ModelViewProtocol) -> None:
"""Set the moved states, and draw the observed state of an emission."""
for agent, value in moves:
setattr(agent, self.state, value)
if self.emission is not None and self.observed is not None:
setattr(agent, self.observed, self._draw_row(self.emission(agent, model)[int(value)], model))
def _move_member(self, agent: AgentViewProtocol, model: ModelViewProtocol) -> None:
"""Move the state of one member, or of the members of its cell, from the states the members before it set."""
if agent.alive:
self._apply(self._get_moves(agent, model), model)
def move(self, model: ModelViewProtocol) -> None:
"""Move the state of every agent of the target that moves in this step, by the schedule.
Args:
model:
The model view of the interaction.
"""
agents = list(model.agents(self._get_role('target')))
if self.schedule == 'synchronous':
moves = [move for agent in agents for move in self._get_moves(agent, model)]
self._apply(moves, model)
return
order = model.rng.permutation(len(agents)) if self.schedule == 'random-sequential' else range(len(agents))
for index in order:
self._move_member(agents[int(index)], model)
def _get_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the jump, the update, the transition or emission law of the state, or the rule of the agent form.
Args:
namespace:
The namespace of the symbols the laws can use, with `namespace.target`.
Returns:
laws:
The law of the move, or the rule that moves the agents.
"""
if self.form == 'agent' and self.schedule == 'member':
return [RuleLaw(self._move_member, phase=self.phase, members=self._get_role('target'))]
if self.form == 'agent':
return [RuleLaw(self.move, phase=self.phase)]
value = getattr(self._get_member(namespace, 'target'), self.state)
if self.matrix is not None:
return self._get_chain(namespace, value)
if self.rate is not None:
rate = self._combine_rates(SR(self.rate(namespace)))
return [self._get_jump(value, self.source, self.destination, rate)]
if self.probability is not None:
probability = self._combine_probabilities(SR(self.probability(namespace)))
moves = self._get_draw() < probability
else:
moves = SR(self.condition(namespace)) if self.condition is not None else SR(0)
if self.after is not None:
moves = both(residence(value) + 1 >= self.after, moves)
return [self._get_step(value, self.source, self.destination, moves)]
Equations