Theta neurons with spike-timing-dependent plasticity
Theta neurons on a random directed network whose synapses learn by spike-timing-dependent plasticity.
- Name:
stdp_network - Field:
biology/neuroscience - Formalism: ordinary differential equations
- References: Song, Miller, and Abbott 2000; Berner et al. 2023
Select this system
from skcomplex.systems import create_simulation, create_system, describe
system = create_system('stdp_network')
entry = describe('stdp_network')
trajectory = create_simulation('stdp_network').prepare().run().trajectory_
Factory
@register_system(
field='biology/neuroscience',
tags=('plasticity', 'stdp', 'adaptive network', 'theta neurons', 'ode'),
references=('Song, Miller, and Abbott 2000', 'Berner et al. 2023'),
defaults={'until': 150.0, 'dt': 0.5},
)
def stdp_network(
n_neurons: int = 20,
*,
p: float = 0.3,
eta: float = 0.5,
spread: float = 0.2,
epsilon: float = 0.05,
A_plus: float = 1.0,
A_minus: float = 1.05,
tau: float = 2.0,
kappa: float = 0.5,
seed: int | None = 0,
) -> System:
"""Return theta neurons on a random directed network whose synapses learn by spike-timing-dependent plasticity.
The neurons fire periodically with excitabilities spread around `eta`, and every synapse is potentiated when its
presynaptic neuron fires shortly before its postsynaptic neuron and depressed when it fires shortly after. Spikes
of neurons that fire at different rates fall at all delays, so the mean change of the efficacies follows the
integral of the window, `(A_plus - A_minus) tau`: a window dominated by depression weakens the synapses and one
dominated by potentiation strengthens them, while the spread of the efficacies grows.
Args:
n_neurons:
The number of neurons.
p:
The probability of every directed synapse.
eta:
The mean excitability.
spread:
The half-width of the uniform spread of the excitabilities.
epsilon:
The learning rate.
A_plus:
The amplitude of potentiation.
A_minus:
The amplitude of depression.
tau:
The time of both windows of the plasticity.
kappa:
The initial efficacy of every synapse.
seed:
The seed of the network, the excitabilities, and the initial phases.
Returns:
system:
The system, with the phases `neurons.theta`, the traces `neurons.a` and `neurons.b`, and the efficacies
`synapses.kappa` on the edges.
Examples:
>>> from skcomplex.systems.biology.neuroscience import stdp_network
>>> stdp_network(10).elements[0].prepare().size_
10
"""
rng = np.random.default_rng(seed)
graph = nx.gnp_random_graph(n_neurons, p, seed=int(rng.integers(2**31)), directed=True)
space = NetworkSpace(graph)
neurons = LatticeElement(
'neurons',
element=TracedThetaNeuronElement,
space=space,
values={
'theta': rng.uniform(0, 2 * np.pi, n_neurons),
'eta': eta + rng.uniform(-spread, spread, n_neurons),
'tau_plus': tau,
'tau_minus': tau,
},
)
synapses = GroupElement('synapses', element=SynapseElement, space=space, on='edges', values={'kappa': kappa})
plasticity = SpikeTimingPlasticityInteraction(
'stdp',
neurons='neurons',
synapses='synapses',
k=max(p * (n_neurons - 1), 1.0),
epsilon=epsilon,
A_plus=A_plus,
A_minus=A_minus,
)
return System(elements=[neurons, synapses], interactions=[plasticity])
Components
Elements:
TracedThetaNeuronElement, fromElementSynapseElement, fromElement
Interactions:
SpikeTimingPlasticityInteraction, fromNeighborInteraction
class TracedThetaNeuronElement(ThetaNeuronElement):
"""Theta neuron that keeps traces of its own spikes, which spike-timing-dependent plasticity reads.
A spike is the smooth pulse `P(theta) = (1 - cos theta)^2 * 2 / 3`, of unit mean over the phase, which peaks as the
phase passes π. The presynaptic trace `a` and the postsynaptic trace `b` rise with every pulse and decay with the
times `tau_plus` and `tau_minus`, `da/dt = P(theta) - a / tau_plus` and `db/dt = P(theta) - b / tau_minus`, so
they hold the recent spikes weighted by the exponential windows of the plasticity.
Args:
label:
The label of the neuron or the group.
visual:
How the element is drawn, or `None` for `default_visual_`.
theta:
The initial phase.
a:
The initial presynaptic trace.
b:
The initial postsynaptic trace.
eta:
The excitability.
tau_plus:
The time of the potentiation window.
tau_minus:
The time of the depression window.
Examples:
>>> from skcomplex.systems.biology.neuroscience import TracedThetaNeuronElement
>>> TracedThetaNeuronElement('neuron').tau_plus
2.0
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
theta: Annotated[float, StateSpec(period=2 * pi)] = 0.0,
a: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
b: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
eta: Annotated[float, ParameterSpec()] = 1.0,
tau_plus: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 2.0,
tau_minus: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 2.0,
) -> None:
super().__init__(label, visual=visual, theta=theta, eta=eta)
self.a = a
self.b = b
self.tau_plus = tau_plus
self.tau_minus = tau_minus
def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return the phase dynamics and the traces of the spikes.
Args:
s:
The namespace of the symbols the laws can use.
Returns:
laws:
The laws of the component.
"""
pulse = 2 * (1 - cos(s.theta)) ** 2 / 3
return [
*super()._get_laws(s),
RateLaw(s.a, pulse - s.a / s.tau_plus),
RateLaw(s.b, pulse - s.b / s.tau_minus),
]
class SynapseElement(Element):
"""Synapse on an edge of a network, with an efficacy that plasticity changes.
Args:
label:
The label of the synapse or of the group on the edges.
visual:
How the element is drawn, or `None` for `default_visual_`.
kappa:
The initial efficacy.
Examples:
>>> from skcomplex.systems.biology.neuroscience import SynapseElement
>>> SynapseElement('synapse', kappa=0.5).kappa
0.5
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
kappa: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.5,
) -> None:
super().__init__(label, visual=visual)
self.kappa = kappa
class SpikeTimingPlasticityInteraction(NeighborInteraction):
"""Pulse coupling along plastic synapses whose efficacy follows the timing of the spikes at both ends.
The presynaptic neuron at the source of an edge drives the postsynaptic neuron at its target with the current
`(1 + cos theta_post) kappa P(theta_pre) / k`, where `k` is the mean in-degree. The efficacy grows when the
postsynaptic neuron fires after the presynaptic one and shrinks when it fires before, `dkappa/dt = epsilon (A_plus
a_pre P(theta_post) - A_minus b_post P(theta_pre))`, the window `W(dt) = A_plus exp(-dt / tau_plus)` for `dt > 0`
and `-A_minus exp(dt / tau_minus)` for `dt < 0` applied to every pair of spikes through the traces. Soft bounds
scale potentiation by `kappa_max - kappa` and depression by `kappa`, so the efficacy stays between 0 and
`kappa_max`.
Args:
label:
The label of the interaction.
neurons:
The path of the group of traced theta neurons on the nodes.
synapses:
The path of the group of synapses on the edges.
k:
The mean in-degree that normalizes the coupling.
epsilon:
The learning rate.
A_plus:
The amplitude of potentiation.
A_minus:
The amplitude of depression.
kappa_max:
The largest efficacy.
Examples:
>>> from skcomplex.systems.biology.neuroscience import SpikeTimingPlasticityInteraction
>>> SpikeTimingPlasticityInteraction('stdp', neurons='neurons', synapses='synapses').A_minus
1.05
"""
def __init__(
self,
label: str,
*,
neurons: Annotated[str, RoleSpec()],
synapses: Annotated[str, RoleSpec()],
k: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
epsilon: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.01,
A_plus: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
A_minus: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.05,
kappa_max: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
super().__init__(label)
self.neurons = neurons
self.synapses = synapses
self.k = k
self.epsilon = epsilon
self.A_plus = A_plus
self.A_minus = A_minus
self.kappa_max = kappa_max
def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return the pulse current into the postsynaptic neuron and the plasticity of the efficacy.
Args:
s:
The namespace of the symbols the laws can use.
Returns:
laws:
The laws of the component.
"""
pre, post, kappa = s.source, s.target, s.edge.kappa
pre_pulse = 2 * (1 - cos(pre.theta)) ** 2 / 3
post_pulse = 2 * (1 - cos(post.theta)) ** 2 / 3
change = s.epsilon * (s.A_plus * pre.a * post_pulse - s.A_minus * post.b * pre_pulse)
bounded = where(change > 0, change * (s.kappa_max - kappa), change * kappa)
return [
RateLaw(post.theta, (1 + cos(post.theta)) * kappa * pre_pulse / s.k),
RateLaw(kappa, bounded / s.kappa_max),
]
Equations
neurons.a[i]: diff(a_neurons(t), t) == 2/3*(cos(theta_neurons) - 1)^2 - a_neurons(t)/tau_plus_neurons
neurons.b[i]: diff(b_neurons(t), t) == 2/3*(cos(theta_neurons) - 1)^2 - b_neurons(t)/tau_minus_neurons