Networks of leaky integrate-and-fire neurons with delays
A random network of leaky integrate-and-fire neurons whose spikes reach their targets after a delay.
- Name:
delayed_lif_network - Field:
biology/neuroscience - Formalism: ordinary differential equations
- References: Brunel and Hakim 1999; arXiv:2112.12147
Factory
@register_model(
field='biology/neuroscience',
tags=('neuron', 'spiking', 'network', 'delay', 'events'),
references=('Brunel and Hakim 1999', 'arXiv:2112.12147'),
defaults={'until': 100.0, 'dt': 0.5},
)
def delayed_lif_network(
n_neurons: int = 10,
*,
probability: float = 0.3,
J: float = 0.5,
D: float = 1.5,
mu: float = 25.0,
random_state: int | np.random.Generator | None = 0,
) -> System:
"""Return a random network of leaky integrate-and-fire neurons whose spikes reach their targets after a delay.
Every neuron is driven above threshold and fires on its own with period `tau ln((mu - reset) / (mu -
threshold))`. Each directed link exists with a probability and moves the target's potential by `J` per spike,
after the delay `D`. Excitation speeds the firing up and inhibition slows it down.
Args:
n_neurons:
The number of neurons.
probability:
The probability of each directed link.
J:
The efficacy of every link in mV.
D:
The transmission delay in ms.
mu:
The constant drive in mV.
random_state:
The seed or generator of the links and the initial potentials.
Returns:
system:
The model, with the potentials `neuron_<i>.v`, the traces `neuron_<i>.s`, and the spikes of neuron `i` as
`spike_<i>` events.
Examples:
>>> from skcomplex.models.biology.neuroscience import delayed_lif_network
>>> len(delayed_lif_network(5).elements)
5
"""
rng = np.random.default_rng(random_state)
potentials = rng.uniform(10.0, 20.0, n_neurons)
neurons = [
SpikingNeuronElement(f'neuron_{index}', v=float(potentials[index]), mu=mu, tau_s=TRACE_TIME)
for index in range(n_neurons)
]
leaks = [
RelaxationInteraction(
f'leak_{index}',
target=neuron.label,
state='v',
goal=lambda s: s.target.mu,
time=lambda s: s.target.tau,
)
for index, neuron in enumerate(neurons)
]
traces = [
RateInteraction(
f'trace_decay_{index}',
target=neuron.label,
state='s',
rate=lambda s: -s.target.s / s.target.tau_s,
)
for index, neuron in enumerate(neurons)
]
spikes = [
SpikeResetInteraction(
f'spike_{index}',
target=neuron.label,
threshold=THRESHOLD,
after=RESET,
kicked='s',
kick=1 / TRACE_TIME,
)
for index, neuron in enumerate(neurons)
]
synapses = [
DelayedSynapseInteraction(
f'synapse_{source}_{target}',
source=f'neuron_{source}',
target=f'neuron_{target}',
J=J,
D=D,
)
for source, target in _random_connections(n_neurons, probability, rng)
]
return System(elements=neurons, interactions=[*leaks, *traces, *spikes, *synapses])
Components
class SpikingNeuronElement(Element):
"""Leaky integrate-and-fire neuron with an exponential synaptic trace of its own spikes.
The potential follows `tau v' = -v + mu` plus synaptic input, and the trace decays as `s' = -s / tau_s`. A spike
reset sets the potential back and raises the trace by `1 / tau_s`, so every spike carries a unit area of trace.
Units are mV and ms. The element holds the states and parameters and has no laws of its own: the leak and trace
interactions of `delayed_lif_network` read its parameters.
Args:
label:
The label of the neuron.
visual:
How the element is drawn, or `None` for `default_visual_`.
v:
The initial potential.
s:
The initial synaptic trace.
tau:
The membrane time constant.
mu:
The constant external drive, the potential the membrane relaxes to.
tau_s:
The decay time of the synaptic trace.
Examples:
>>> from skcomplex.models.biology.neuroscience import SpikingNeuronElement
>>> SpikingNeuronElement('neuron').tau_s
2.0
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
v: Annotated[float, StateSpec()] = 0.0,
s: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
tau: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 20.0,
mu: Annotated[float, ParameterSpec()] = 25.0,
tau_s: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 2.0,
) -> None:
super().__init__(label, visual=visual)
self.v = v
self.s = s
self.tau = tau
self.mu = mu
self.tau_s = tau_s
class RelaxationInteraction(RateInteraction):
"""Relaxation `(goal(ns) - x) / time(ns)` of a state `x` of a target towards a goal, with optional noise.
With a constant goal and time and an additive noise of amplitude `sigma`, the state follows the Ornstein-Uhlenbeck
process `dx = (goal - x) / time dt + sigma dW`, whose stationary variance is `sigma**2 time / 2`. Without noise it
is first-order smoothing, or adaptive expectation, towards the goal. A goal or a time function takes the namespace
of the interaction, as for `RateInteraction`, and a number stands for a constant. A goal that is a number applies
to every component of a vector state.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose state relaxes.
state:
The name of the state of the target that relaxes.
goal:
A function of the namespace that returns the goal, or a number.
time:
A function of the namespace that returns the relaxation time, or a positive number.
noise:
A function of the namespace that returns the amplitude of the noise, a number, or `None` for no noise.
Raises:
ValidationError:
The state is not a name, or a function is neither callable nor a number.
Examples:
>>> from skcomplex.core.elements import CompartmentElement
>>> from skcomplex.core.interactions import RelaxationInteraction
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> smoothing = RelaxationInteraction('smoothing', target='belief', state='n', goal=10.0, time=0.5)
>>> system = System(elements=[CompartmentElement('belief', n=0.0)], interactions=[smoothing])
>>> trajectory = Simulation(system, dt=0.01).check().prepare().run(until=5.0).trajectory_
>>> round(float(trajectory.values('belief.n')[-1]), 3)
10.0
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
state: str = 'x',
goal: Callable[[SimpleNamespace], object] | float = 0.0,
time: Callable[[SimpleNamespace], object] | float = 1.0,
noise: Callable[[SimpleNamespace], object] | float | None = None,
) -> None:
super().__init__(label, target=target, state=state, noise=noise)
self.goal = goal
self.time = time
def _check_functions(self) -> None:
"""Raise an error for a state that is not a name, or a function that is neither callable nor a number."""
super()._check_functions()
for name in ('goal', 'time'):
value = getattr(self, name)
if not (callable(value) or isinstance(value, int | float)):
error_msg = f'The {name} of {self.label!r} is {value!r}. Give a function of the namespace or a number.'
raise ValidationError(error_msg)
def _get_goal(self, namespace: SimpleNamespace) -> object:
"""Return the goal, the hook a subclass overrides to derive it from its parameters."""
return self._evaluate(self.goal, namespace)
def _get_time(self, namespace: SimpleNamespace) -> object:
"""Return the relaxation time, the hook a subclass overrides to derive it from its parameters."""
return self._evaluate(self.time, namespace)
def _get_rate(self, namespace: SimpleNamespace) -> object:
"""Return the relaxation `(goal - x) / time` of the state.
Args:
namespace:
The namespace of the symbols the rate can use.
Returns:
rate:
The rate of the state.
"""
state = getattr(self._get_member(namespace, 'target'), self.state)
goal = self._broadcast(self._get_goal(namespace), state)
return (goal - state) / self._get_time(namespace)
class RateInteraction(Interaction):
"""Rate `rate(ns)` and noise `noise(ns) dW` added to the highest time derivative of a state of a target.
On a state of order 1 the rate is its first derivative, and on a state of order `k` it adds to the `k`-th
derivative, so rates on one state superpose. With a divisor `d` the rate and the noise are divided by `d`, the name
of a state or parameter of the target, such as an inertia, a friction, or a capacitance, which backends can read.
A rate, a noise, or a divisor function takes the namespace of the interaction, with the namespaces `target` and
`sources` of the bound elements, the time `t`, and the symbols of the parameters a subclass declares. A number
stands for a constant, and a number on a vector state applies to every component. The noise is a Wiener process per
component and member by default, a process of a `WienerProcessSpec` to share or correlate it across states and
members, or a fractional Brownian motion of a Hurst index.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose state the rate changes.
state:
The name of the state of the target that the rate changes.
rate:
A function of the namespace that returns the rate, or a number.
sources:
The paths of the elements or groups the rate and the noise read but do not change.
noise:
A function of the namespace that returns the amplitude of the noise, a number, or `None` for no noise.
divisor:
The name of the state or parameter of the target that divides the rate and the noise, a function of the
namespace of the target that returns the divisor, such as a moment of inertia `m l^2`, or `None`.
process:
The Wiener processes that drive the noise, or `None` for one independent process per component and member.
index:
The process of `process` that drives the noise, or `None` for every process, one per component of a vector
state.
convention:
The stochastic calculus of the noise, `'ito'`, `'stratonovich'`, or `'anti-ito'`.
hurst:
The Hurst index in `(0, 1)` of a fractional Brownian motion that drives the noise, a number or a function
of the namespace that returns a parameter, or `None` for a Wiener process.
conserved:
Whether the noise is the divergence of a random flux on the links of the lattice of a field, so it keeps the
total of the state.
Raises:
ValidationError:
The state is not a name, a function is neither callable nor a number, the convention is unknown, the Hurst
index is outside `(0, 1)`, or the noise options cannot combine.
Examples:
>>> from skcomplex.core.elements import CompartmentElement
>>> from skcomplex.core.interactions import RateInteraction
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> decay = RateInteraction('decay', target='stock', state='n', rate=lambda s: -0.5 * s.target.n)
>>> system = System(elements=[CompartmentElement('stock', n=8.0)], interactions=[decay])
>>> trajectory = Simulation(system, dt=0.01).check().prepare().run(until=2.0).trajectory_
>>> round(float(trajectory.values('stock.n')[-1]), 3)
2.943
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
state: str = 'x',
rate: Callable[[SimpleNamespace], object] | float = 0.0,
sources: Annotated[Sequence[str], RoleSpec(variadic=True)] = (),
noise: Callable[[SimpleNamespace], object] | float | None = None,
divisor: str | Callable[[SimpleNamespace], object] | None = None,
process: Annotated[WienerProcessSpec | None, RoleSpec(noise=True)] = None,
index: int | None = None,
convention: str = 'ito',
hurst: Callable[[SimpleNamespace], object] | float | None = None,
conserved: bool = False,
) -> None:
super().__init__(label)
self.target = target
self.state = state
self.rate = rate
self.sources = sources
self.noise = noise
self.divisor = divisor
self.process = process
self.index = index
self.convention = convention
self.hurst = hurst
self.conserved = conserved
def _check_functions(self) -> None:
"""Raise an error for a state that is not a name, or a function that is neither callable nor a number."""
if not isinstance(self.state, str) or not self.state.isidentifier():
error_msg = (
f'{type(self).__name__} {self.label!r} changes the state {self.state!r}, which is not a name of a '
'state.'
)
raise ValidationError(error_msg)
for name in ('rate', 'noise'):
value = getattr(self, name)
if not (value is None or callable(value) or isinstance(value, int | float)):
error_msg = f'The {name} of {self.label!r} is {value!r}. Give a function of the namespace or a number.'
raise ValidationError(error_msg)
if not (
self.divisor is None
or callable(self.divisor)
or (isinstance(self.divisor, str) and self.divisor.isidentifier())
):
error_msg = (
f'The divisor of {self.label!r} is {self.divisor!r}. Give the name of a state or parameter, or a '
'function of the namespace of the target.'
)
raise ValidationError(error_msg)
def _check_noise(self) -> None:
"""Raise an error for an unknown convention, a Hurst index outside `(0, 1)`, or noise options that conflict."""
if self.convention not in CONVENTIONS:
error_msg = f'Unknown noise convention {self.convention!r}. Use one of {", ".join(CONVENTIONS)}.'
raise ValidationError(error_msg)
if not (
self.hurst is None or callable(self.hurst) or (isinstance(self.hurst, int | float) and 0 < self.hurst < 1)
):
error_msg = f'The Hurst index of {self.label!r} is {self.hurst!r}. Give a number in (0, 1) or a function.'
raise ValidationError(error_msg)
if self.index is not None and self.process is None:
error_msg = f'{self.label!r} picks the process {self.index!r} but binds no processes. Give `process`.'
raise ValidationError(error_msg)
options = [
name
for name, given in (('hurst', self.hurst), ('process', self.process), ('conserved', self.conserved))
if given
]
if len(options) > 1:
error_msg = f'The noise of {self.label!r} takes {" and ".join(options)}. Give at most one of them.'
raise ValidationError(error_msg)
def _check(self) -> None:
"""Raise an error for a state, functions, or noise options the interaction cannot take.
Raises:
ValidationError:
The type declares states, the state is not a name, a function is neither callable nor a number, the
convention is unknown, the Hurst index is outside `(0, 1)`, or the noise options cannot combine.
"""
super()._check()
self._check_functions()
self._check_noise()
def _evaluate(self, function: Callable[[SimpleNamespace], object] | float, namespace: SimpleNamespace) -> object:
"""Return the value of a function of the namespace, or a number as it is."""
return function(namespace) if callable(function) else function
def _broadcast(self, value: object, state: object) -> object:
"""Return a value with one component per component of a vector state, repeating a number."""
if isinstance(state, FreeModuleElement) and not isinstance(value, FreeModuleElement | list | tuple):
return vector(SR, [value] * len(state))
if isinstance(value, list | tuple):
return vector(SR, value)
return value
def _divide(self, value: object, namespace: SimpleNamespace) -> object:
"""Return a value divided by the divisor of the target, or the value without a divisor."""
if self.divisor is None:
return value
divisor = (
self.divisor(self._get_member(namespace, 'target'))
if callable(self.divisor)
else getattr(self._get_member(namespace, 'target'), self.divisor)
)
return operator.truediv(value, divisor)
def _get_rate(self, namespace: SimpleNamespace) -> object:
"""Return the rate before the divisor, the hook a subclass overrides to derive it from its parameters."""
return self._evaluate(self.rate, namespace)
def _get_noise(self, namespace: SimpleNamespace) -> object:
"""Return the noise amplitude before the divisor, the hook a subclass overrides, or `None` for no noise."""
return None if self.noise is None else self._evaluate(self.noise, namespace)
def _get_noise_law(self, state: object, amplitude: object, namespace: SimpleNamespace) -> BasePart:
"""Return the noise law of an amplitude, by the process, the convention, and the Hurst index."""
if self.hurst is not None:
hurst = self._evaluate(self.hurst, namespace)
return FractionalNoiseLaw(state, amplitude, hurst=hurst, convention=self.convention)
process = None
if self.process is not None:
handle = self._get_member(namespace, 'process').dW
process = handle if self.index is None else handle[self.index]
return NoiseLaw(state, amplitude, convention=self.convention, process=process, conserved=self.conserved)
def _get_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the rate law of the state and its noise law, both divided by the divisor.
Args:
namespace:
The namespace of the symbols the laws can use, with the namespaces of the target and the sources.
Returns:
laws:
The rate law, and the noise law when the interaction has noise.
"""
state = getattr(self._get_member(namespace, 'target'), self.state)
rate = self._divide(self._broadcast(self._get_rate(namespace), state), namespace)
laws: list[BasePart] = [RateLaw(state, rate)]
noise = self._get_noise(namespace)
if noise is not None:
laws.append(self._get_noise_law(state, self._divide(noise, namespace), namespace))
return laws
class SpikeResetInteraction(CrossingInteraction):
"""Spike of a neuron, which resets its potential when it crosses a threshold upwards and kicks a second state.
When the potential `v` rises through `threshold`, it is set to `after`, and the state `kicked`, such as a recovery
variable or a synaptic trace, grows by `kick`. Every spike is an event labelled by the interaction.
Args:
label:
The label of the interaction and of its events.
target:
The path of the neuron.
state:
The name of the potential.
threshold:
The potential of a spike.
after:
The potential after a spike.
kicked:
The name of the state a spike increments, or `None` for none.
kick:
The increment of the kicked state at a spike.
Raises:
ValidationError:
The direction is not -1, 0, or 1, or the guard or a new value is not a function.
Examples:
>>> from skcomplex.core.elements import CompartmentElement
>>> from skcomplex.core.interactions import RateInteraction
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.models.biology.neuroscience import SpikeResetInteraction
>>> drive = RateInteraction('drive', target='neuron', state='n', rate=1.0)
>>> spike = SpikeResetInteraction('spike', target='neuron', state='n', threshold=1.0, after=0.0)
>>> system = System(elements=[CompartmentElement('neuron')], interactions=[drive, spike])
>>> trajectory = Simulation(system, dt=0.25).check().prepare().run(until=2.1).trajectory_
>>> [round(event.time, 6) for event in trajectory.events]
[1.0, 2.0]
"""
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec()],
state: str = 'v',
threshold: Annotated[float, ParameterSpec()] = 1.0,
after: Annotated[float, ParameterSpec()] = 0.0,
kicked: str | None = None,
kick: Annotated[float, ParameterSpec()] = 0.0,
) -> None:
super().__init__(
label,
target=target,
guard=self._get_guard,
reset=self._get_resets(state, kicked),
direction=1,
)
self.state = state
self.threshold = threshold
self.after = after
self.kicked = kicked
self.kick = kick
def _get_guard(self, s: SimpleNamespace) -> Expression:
"""Return the potential above the threshold."""
return getattr(s.target, self.state) - s.threshold
def _get_resets(self, state: str, kicked: str | None) -> dict[str, Callable[[SimpleNamespace], object]]:
"""Return the new value of the potential and of the kicked state after a spike."""
resets: dict[str, Callable[[SimpleNamespace], object]] = {state: lambda s: s.after}
if kicked is not None:
resets[kicked] = lambda s: getattr(s.target, kicked) + s.kick
return resets
class DelayedSynapseInteraction(RateInteraction):
"""Synapse that adds the presynaptic trace from a transmission delay ago to the postsynaptic potential.
The target gains `J s_source(t - D)`, so every presynaptic spike moves the postsynaptic potential by `J` in
total, spread over the decay time of the trace and arriving after the delay `D`.
Args:
label:
The label of the interaction.
source:
The path of the presynaptic neuron.
target:
The path of the postsynaptic neuron.
J:
The efficacy in mV, positive for excitation and negative for inhibition.
D:
The transmission delay in ms.
Examples:
>>> from skcomplex.models.biology.neuroscience import DelayedSynapseInteraction
>>> DelayedSynapseInteraction('synapse', source='a', target='b', J=-1.0).D
1.5
"""
def __init__(
self,
label: str,
*,
source: Annotated[str, RoleSpec()],
target: Annotated[str, RoleSpec()],
J: Annotated[float, ParameterSpec()] = 0.5,
D: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.5,
) -> None:
super().__init__(label, target=target, state='v')
self.source = source
self.J = J
self.D = D
def _get_rate(self, namespace: SimpleNamespace) -> Expression:
"""Return the delayed synaptic input `J s_source(t - D)`.
Args:
namespace:
The namespace of the symbols the rate can use, with the namespaces of the source and the target.
Returns:
rate:
The rate of the postsynaptic potential.
"""
return namespace.J * delayed(namespace.source.s, namespace.D)
Equations
neuron_0.v: diff(v_neuron_0(t), t) == J_synapse_4_0*s_neuron_4(-D_synapse_4_0 + t) + J_synapse_5_0*s_neuron_5(-D_synapse_5_0 + t) + J_synapse_6_0*s_neuron_6(-D_synapse_6_0 + t) + (mu_neuron_0 - v_neuron_0(t))/tau_neuron_0
neuron_0.s: diff(s_neuron_0(t), t) == -s_neuron_0(t)/tau_s_neuron_0
neuron_1.v: diff(v_neuron_1(t), t) == J_synapse_9_1*s_neuron_9(-D_synapse_9_1 + t) + (mu_neuron_1 - v_neuron_1(t))/tau_neuron_1
neuron_1.s: diff(s_neuron_1(t), t) == -s_neuron_1(t)/tau_s_neuron_1
neuron_2.v: diff(v_neuron_2(t), t) == J_synapse_0_2*s_neuron_0(-D_synapse_0_2 + t) + J_synapse_1_2*s_neuron_1(-D_synapse_1_2 + t) + J_synapse_4_2*s_neuron_4(-D_synapse_4_2 + t) + J_synapse_6_2*s_neuron_6(-D_synapse_6_2 + t) + (mu_neuron_2 - v_neuron_2(t))/tau_neuron_2
neuron_2.s: diff(s_neuron_2(t), t) == -s_neuron_2(t)/tau_s_neuron_2
neuron_3.v: diff(v_neuron_3(t), t) == J_synapse_1_3*s_neuron_1(-D_synapse_1_3 + t) + J_synapse_5_3*s_neuron_5(-D_synapse_5_3 + t) + (mu_neuron_3 - v_neuron_3(t))/tau_neuron_3
neuron_3.s: diff(s_neuron_3(t), t) == -s_neuron_3(t)/tau_s_neuron_3
neuron_4.v: diff(v_neuron_4(t), t) == J_synapse_0_4*s_neuron_0(-D_synapse_0_4 + t) + J_synapse_5_4*s_neuron_5(-D_synapse_5_4 + t) + (mu_neuron_4 - v_neuron_4(t))/tau_neuron_4
neuron_4.s: diff(s_neuron_4(t), t) == -s_neuron_4(t)/tau_s_neuron_4
neuron_5.v: diff(v_neuron_5(t), t) == J_synapse_2_5*s_neuron_2(-D_synapse_2_5 + t) + J_synapse_6_5*s_neuron_6(-D_synapse_6_5 + t) + J_synapse_9_5*s_neuron_9(-D_synapse_9_5 + t) + (mu_neuron_5 - v_neuron_5(t))/tau_neuron_5
neuron_5.s: diff(s_neuron_5(t), t) == -s_neuron_5(t)/tau_s_neuron_5
neuron_6.v: diff(v_neuron_6(t), t) == J_synapse_0_6*s_neuron_0(-D_synapse_0_6 + t) + J_synapse_4_6*s_neuron_4(-D_synapse_4_6 + t) + J_synapse_8_6*s_neuron_8(-D_synapse_8_6 + t) + (mu_neuron_6 - v_neuron_6(t))/tau_neuron_6
neuron_6.s: diff(s_neuron_6(t), t) == -s_neuron_6(t)/tau_s_neuron_6
neuron_7.v: diff(v_neuron_7(t), t) == J_synapse_5_7*s_neuron_5(-D_synapse_5_7 + t) + (mu_neuron_7 - v_neuron_7(t))/tau_neuron_7
neuron_7.s: diff(s_neuron_7(t), t) == -s_neuron_7(t)/tau_s_neuron_7
neuron_8.v: diff(v_neuron_8(t), t) == J_synapse_4_8*s_neuron_4(-D_synapse_4_8 + t) + J_synapse_5_8*s_neuron_5(-D_synapse_5_8 + t) + (mu_neuron_8 - v_neuron_8(t))/tau_neuron_8
neuron_8.s: diff(s_neuron_8(t), t) == -s_neuron_8(t)/tau_s_neuron_8
neuron_9.v: diff(v_neuron_9(t), t) == J_synapse_0_9*s_neuron_0(-D_synapse_0_9 + t) + (mu_neuron_9 - v_neuron_9(t))/tau_neuron_9
neuron_9.s: diff(s_neuron_9(t), t) == -s_neuron_9(t)/tau_s_neuron_9