Neuroscience

Neuroscience: spiking and conductance-based neurons, populations, networks, excitable networks, and synapses.

AdaptrodeElement(label, visual=None, *, w0=0.0, w1=0.0, w2=0.0, w_max=1.0, alpha0=0.25, alpha1=0.05, alpha2=0.01, delta0=0.1, delta1=0.02, delta2=0.002, amplitude=1.0, period=10.0, width=3.0, stop=150.0)

Bases: Element

Adaptrode of Mobus, a synaptic weight with three levels of memory that learn on increasing time scales.

At every step the fast weight w0 rises toward w_max in proportion to the input x and decays toward the slower weight w1, which rises toward w0 and decays toward the slowest weight w2, which rises toward w1 and decays toward zero. The input is a train of pulses of height amplitude, width steps long, every period steps, until the step stop, so the weights record recent, medium-term, and long-term activity. The element holds the states and parameters and has no laws of its own: the maps of adaptrode, one per level, read its parameters.

Parameters:

Name	Type	Description	Default
label	str	The label of the synapse.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
w0	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial fast weight.	0.0
w1	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial medium weight.	0.0
w2	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial slow weight.	0.0
w_max	Annotated[float, ParameterSpec(domain=open(0, oo))]	The largest weight.	1.0
alpha0	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The learning rate of the fast weight.	0.25
alpha1	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The learning rate of the medium weight.	0.05
alpha2	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The learning rate of the slow weight.	0.01
delta0	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The decay rate of the fast weight.	0.1
delta1	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The decay rate of the medium weight.	0.02
delta2	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The decay rate of the slow weight.	0.002
amplitude	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The height of the input pulses.	1.0
period	Annotated[float, ParameterSpec(domain=open(0, oo))]	The number of steps between the starts of two pulses.	10.0
width	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The number of steps of each pulse.	3.0
stop	Annotated[float, ParameterSpec()]	The step after which the input stays off.	150.0

Examples:

>>> from skcomplex.systems.biology.neuroscience import AdaptrodeElement
>>> AdaptrodeElement('synapse').w_max
1.0

Source code in src/skcomplex/systems/biology/neuroscience/_synapses.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    w0: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    w1: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    w2: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    w_max: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    alpha0: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.25,
    alpha1: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.05,
    alpha2: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.01,
    delta0: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
    delta1: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.02,
    delta2: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.002,
    amplitude: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    period: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 10.0,
    width: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 3.0,
    stop: Annotated[float, ParameterSpec()] = 150.0,
) -> None:
    super().__init__(label, visual=visual)
    self.w0 = w0
    self.w1 = w1
    self.w2 = w2
    self.w_max = w_max
    self.alpha0 = alpha0
    self.alpha1 = alpha1
    self.alpha2 = alpha2
    self.delta0 = delta0
    self.delta1 = delta1
    self.delta2 = delta2
    self.amplitude = amplitude
    self.period = period
    self.width = width
    self.stop = stop

AssociativeSiteElement(label, visual=None, *, x=0.0, phi=1.0, x_c=0.85, gamma_plus=0.005, gamma_minus=0.005)

Bases: Element

Site of a dense homogeneous associative network, an activity coupled to a slow reservoir.

The reservoir phi of an active site, x > x_c, depletes as dphi/dt = -gamma_minus phi, and the reservoir of an inactive site recovers as dphi/dt = gamma_plus (1 - phi) (1 - x / x_c). The network couples the activities.

Parameters:

Name	Type	Description	Default
label	str	The label of the site or of the group of sites.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec(domain=UNIT)]	The initial activity.	0.0
phi	Annotated[float, StateSpec(domain=UNIT)]	The initial reservoir.	1.0
x_c	Annotated[float, ParameterSpec(domain=open(0, 1))]	The activity above which a site counts as active.	0.85
gamma_plus	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The recovery rate of the reservoir.	0.005
gamma_minus	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The depletion rate of the reservoir.	0.005

Examples:

>>> from skcomplex.systems.biology.neuroscience import AssociativeSiteElement
>>> AssociativeSiteElement('sites').x_c
0.85

Source code in src/skcomplex/systems/biology/neuroscience/_associative.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec(domain=UNIT)] = 0.0,
    phi: Annotated[float, StateSpec(domain=UNIT)] = 1.0,
    x_c: Annotated[float, ParameterSpec(domain=RealSet.open(0, 1))] = 0.85,
    gamma_plus: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.005,
    gamma_minus: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.005,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.phi = phi
    self.x_c = x_c
    self.gamma_plus = gamma_plus
    self.gamma_minus = gamma_minus

AxonElement(label, visual=None, *, tau=1.0, tau0=1.0)

Bases: Element

Axon on an edge of a network, whose conduction delay myelination changes.

Parameters:

Name	Type	Description	Default
label	str	The label of the axon or of the group on the edges.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
tau	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial conduction delay.	1.0
tau0	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The baseline delay that the delay relaxes to.	1.0

Examples:

>>> from skcomplex.systems.biology.neuroscience import AxonElement
>>> AxonElement('axon', tau0=2.0).tau0
2.0

Source code in src/skcomplex/systems/biology/neuroscience/_plastic_networks.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    tau: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    tau0: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.tau = tau
    self.tau0 = tau0

CliqueCompetitionInteraction(label, *, sites, w=((0.0,),), z=1.0, phi_w=0.7, phi_z=0.15, width=0.05, f_w_min=0.1, f_z_min=0.0)

Bases: Interaction

Competition of the cliques of a winners-take-all network, gated by the reservoirs of the sites.

Site i grows by dx_i/dt = (1 - x_i) r_i for a positive input r_i and decays by dx_i/dt = x_i r_i for a negative one, with r_i = sum_j [f_w(phi_i) w_ij + z_ij f_z(phi_j)] x_j. The excitatory links w_ij > 0 join the sites of every clique, and every other pair inhibits with z_ij = -|z|. The reservoir functions f(phi) = f_min + (1 - f_min) (atan((phi - phi_c) / width) - atan(-phi_c / width)) / (atan((1 - phi_c) / width) - atan(-phi_c / width)) switch off the excitation of a site whose own reservoir is low and the inhibition by a site whose reservoir is low.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
sites	Annotated[str, RoleSpec()]	The path of the group of sites.	required
w	Sequence[Sequence[float]]	The excitatory links, a symmetric matrix with zeros where sites inhibit each other.	((0.0,),)
z	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The strength of the inhibition.	1.0
phi_w	Annotated[float, ParameterSpec(domain=open(0, 1))]	The inflection point of the reservoir function of the excitation.	0.7
phi_z	Annotated[float, ParameterSpec(domain=open(0, 1))]	The inflection point of the reservoir function of the inhibition.	0.15
width	Annotated[float, ParameterSpec(domain=open(0, oo))]	The width of the reservoir functions.	0.05
f_w_min	Annotated[float, ParameterSpec(domain=closed(0, 1))]	The floor of the reservoir function of the excitation.	0.1
f_z_min	Annotated[float, ParameterSpec(domain=closed(0, 1))]	The floor of the reservoir function of the inhibition.	0.0

Examples:

>>> from skcomplex.systems.biology.neuroscience import CliqueCompetitionInteraction
>>> CliqueCompetitionInteraction('competition', sites='sites', w=((0.0, 0.1), (0.1, 0.0))).z
1.0

Source code in src/skcomplex/systems/biology/neuroscience/_associative.py

def __init__(
    self,
    label: str,
    *,
    sites: Annotated[str, RoleSpec()],
    w: Sequence[Sequence[float]] = ((0.0,),),
    z: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    phi_w: Annotated[float, ParameterSpec(domain=RealSet.open(0, 1))] = 0.7,
    phi_z: Annotated[float, ParameterSpec(domain=RealSet.open(0, 1))] = 0.15,
    width: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.05,
    f_w_min: Annotated[float, ParameterSpec(domain=RealSet.closed(0, 1))] = 0.1,
    f_z_min: Annotated[float, ParameterSpec(domain=RealSet.closed(0, 1))] = 0.0,
) -> None:
    super().__init__(label)
    self.sites = sites
    self.w = w
    self.z = z
    self.phi_w = phi_w
    self.phi_z = phi_z
    self.width = width
    self.f_w_min = f_w_min
    self.f_z_min = f_z_min

DelayAdaptationInteraction(label, *, source, target, axon, kappa=0.1, epsilon=0.01, rho=1.0, width=0.05, shortest=0.1)

Bases: Interaction

Phase coupling through conduction delays that adapt to the phase differences of the units they join.

The target unit reads the phase of the source unit along the axon as it was one conduction delay ago, dtheta_i/dt += kappa sin(theta_j(t - tau_ij) - theta_i). The delay relaxes to its baseline and is pushed by the phase difference, dtau_ij/dt = epsilon H(tau_ij) (-(tau_ij - tau0_ij) + rho sin(theta_j - theta_i)), where the smooth step H(tau) = (1 + tanh(tau / width)) / 2 keeps the delay nonnegative. The coupling reads the delay no shorter than shortest, which bounds the steps of the solver.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
source	Annotated[str, RoleSpec()]	The path of the unit that sends the signal.	required
target	Annotated[str, RoleSpec()]	The path of the unit that receives the signal.	required
axon	Annotated[str, RoleSpec()]	The path of the axon that conducts it.	required
kappa	Annotated[float, ParameterSpec()]	The coupling strength of an edge.	0.1
epsilon	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The adaptation rate of the delays.	0.01
rho	Annotated[float, ParameterSpec()]	The strength of the push of the phase difference on the delay.	1.0
width	Annotated[float, ParameterSpec(domain=open(0, oo))]	The width of the smooth step.	0.05
shortest	float	The shortest delay the coupling reads.	0.1

Examples:

>>> from skcomplex.systems.biology.neuroscience import DelayAdaptationInteraction
>>> DelayAdaptationInteraction('myelination', source='a', target='b', axon='ab', rho=0.5).rho
0.5

Source code in src/skcomplex/systems/biology/neuroscience/_plastic_networks.py

def __init__(
    self,
    label: str,
    *,
    source: Annotated[str, RoleSpec()],
    target: Annotated[str, RoleSpec()],
    axon: Annotated[str, RoleSpec()],
    kappa: Annotated[float, ParameterSpec()] = 0.1,
    epsilon: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.01,
    rho: Annotated[float, ParameterSpec()] = 1.0,
    width: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.05,
    shortest: float = 0.1,
) -> None:
    super().__init__(label)
    self.source = source
    self.target = target
    self.axon = axon
    self.kappa = kappa
    self.epsilon = epsilon
    self.rho = rho
    self.width = width
    self.shortest = shortest

DelayedSynapseInteraction(label, *, source, target, J=0.5, D=1.5)

Bases: 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.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
source	Annotated[str, RoleSpec()]	The path of the presynaptic neuron.	required
target	Annotated[str, RoleSpec()]	The path of the postsynaptic neuron.	required
J	Annotated[float, ParameterSpec()]	The efficacy in mV, positive for excitation and negative for inhibition.	0.5
D	Annotated[float, ParameterSpec(domain=open(0, oo))]	The transmission delay in ms.	1.5

Examples:

>>> from skcomplex.systems.biology.neuroscience import DelayedSynapseInteraction
>>> DelayedSynapseInteraction('synapse', source='a', target='b', J=-1.0).D
1.5

Source code in src/skcomplex/systems/biology/neuroscience/_networks.py

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

ExcitableElement(label, visual=None, *, phase=0, p=0.1, stimulus=0.0)

Bases: 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.

Parameters:

Name	Type	Description	Default
label	str	The label of the element or the site group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
phase	Annotated[float, StateSpec(domain=NN)]	The initial state, 0 for quiescent, 1 for excited, and 2 to n_states - 1 for refractory.	0
p	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The probability that an excited neighbour excites the element in one step.	0.1
stimulus	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The rate r of the external stimulus per step.	0.0

Examples:

>>> from skcomplex.systems.biology.neuroscience import ExcitableElement
>>> ExcitableElement('neuron').p
0.1

Source code in src/skcomplex/systems/biology/neuroscience/_excitable.py

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

FitzHughNagumoNeuronElement(label, visual=None, *, v=-1.2, w=-0.6)

Bases: Element

FitzHugh-Nagumo neuron, with a fast potential v and a slow recovery w.

The neuron has no laws of its own. FitzHugh-Nagumo kinetics change it.

Parameters:

Name	Type	Description	Default
label	str	The label of the neuron.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
v	Annotated[float, StateSpec()]	The initial potential.	-1.2
w	Annotated[float, StateSpec()]	The initial recovery.	-0.6

Examples:

>>> from skcomplex.systems.biology.neuroscience import FitzHughNagumoNeuronElement
>>> FitzHughNagumoNeuronElement('neuron').w
-0.6

Source code in src/skcomplex/systems/biology/neuroscience/_conductances.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    v: Annotated[float, StateSpec()] = -1.2,
    w: Annotated[float, StateSpec()] = -0.6,
) -> None:
    super().__init__(label, visual=visual)
    self.v = v
    self.w = w

GlobalInhibitionInteraction(label, *, oscillators, inhibitor, W_z=1.5, theta_z=0.1)

Bases: CouplingInteraction

Inhibition of every oscillator of a lattice by a global inhibitor, once the inhibitor exceeds its threshold.

Every oscillator gets -W_z H(z - theta_z) from the inhibitor, a smooth step of width kappa of the oscillator. It is the coupling of the members of a group to a hub element, and InhibitorDriveInteraction drives the hub.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
oscillators	Annotated[str, RoleSpec()]	The path of the lattice of oscillators.	required
inhibitor	Annotated[str, RoleSpec()]	The path of the inhibitor.	required
W_z	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The strength of the inhibition.	1.5
theta_z	Annotated[float, ParameterSpec()]	The threshold of the inhibitor.	0.1

Examples:

>>> from skcomplex.systems.biology.neuroscience import GlobalInhibitionInteraction
>>> GlobalInhibitionInteraction('inhibition', inhibitor='z', oscillators='grid').W_z
1.5

Source code in src/skcomplex/systems/biology/neuroscience/_networks.py

def __init__(
    self,
    label: str,
    *,
    oscillators: Annotated[str, RoleSpec()],
    inhibitor: Annotated[str, RoleSpec()],
    W_z: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.5,
    theta_z: Annotated[float, ParameterSpec()] = 0.1,
) -> None:
    super().__init__(
        label,
        target=oscillators,
        source=inhibitor,
        write='x',
        function=self._get_inhibition,
        mode='group',
    )
    self.oscillators = oscillators
    self.inhibitor = inhibitor
    self.W_z = W_z
    self.theta_z = theta_z

HodgkinHuxleyElement(label, visual=None, *, v=-65.0, n=0.3177, m=0.0529, h=0.5961, C=1.0, g_na=120.0, g_k=36.0, g_l=0.3, e_na=50.0, e_k=-77.0, e_l=-54.387, I=10.0)

Bases: Element

Squid giant axon membrane of Hodgkin and Huxley, with sodium, potassium, and leak currents.

The potential follows C v' = I - g_Na m³ h (v - E_Na) - g_K n⁴ (v - E_K) - g_L (v - E_L), and each gate x of n, m, and h relaxes as x' = alpha_x(v) (1 - x) - beta_x(v) x, in the modern convention with rest near -65 mV. Upward crossings of 0 mV are recorded as spike events. Units are mV, ms, mS/cm², and µA/cm². The element holds the states and parameters and has no laws of its own: the current, gating, and crossing interactions of hodgkin_huxley read its parameters.

Parameters:

Name	Type	Description	Default
label	str	The label of the neuron.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
v	Annotated[float, StateSpec()]	The initial membrane potential.	-65.0
n	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial potassium activation.	0.3177
m	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial sodium activation.	0.0529
h	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial sodium inactivation.	0.5961
C	Annotated[float, ParameterSpec(domain=open(0, oo))]	The membrane capacitance.	1.0
g_na	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The maximal sodium conductance.	120.0
g_k	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The maximal potassium conductance.	36.0
g_l	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The leak conductance.	0.3
e_na	Annotated[float, ParameterSpec()]	The sodium reversal potential.	50.0
e_k	Annotated[float, ParameterSpec()]	The potassium reversal potential.	-77.0
e_l	Annotated[float, ParameterSpec()]	The leak reversal potential.	-54.387
I	Annotated[float, ParameterSpec()]	The applied current.	10.0

Examples:

>>> from skcomplex.systems.biology.neuroscience import HodgkinHuxleyElement
>>> HodgkinHuxleyElement('axon', I=10.0).g_na
120.0

Source code in src/skcomplex/systems/biology/neuroscience/_conductances.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    v: Annotated[float, StateSpec()] = -65.0,
    n: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.3177,
    m: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0529,
    h: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.5961,
    C: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    g_na: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 120.0,
    g_k: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 36.0,
    g_l: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.3,
    e_na: Annotated[float, ParameterSpec()] = 50.0,
    e_k: Annotated[float, ParameterSpec()] = -77.0,
    e_l: Annotated[float, ParameterSpec()] = -54.387,
    I: Annotated[float, ParameterSpec()] = 10.0,
) -> None:
    super().__init__(label, visual=visual)
    self.v = v
    self.n = n
    self.m = m
    self.h = h
    self.C = C
    self.g_na = g_na
    self.g_k = g_k
    self.g_l = g_l
    self.e_na = e_na
    self.e_k = e_k
    self.e_l = e_l
    self.I = I

InhibitorDriveInteraction(label, *, inhibitor, oscillators)

Bases: CouplingInteraction

Drive of a global inhibitor by the largest activity m H(x - theta_x) of the oscillators of a lattice.

The inhibitor gets phi max_j m_j H(x_j - theta_x), with the rate phi of the inhibitor and a smooth step of width kappa of every oscillator, so it is driven when at least one stimulated oscillator is active. It is the coupling of a hub element to a group by the largest value over the group.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
inhibitor	Annotated[str, RoleSpec()]	The path of the inhibitor.	required
oscillators	Annotated[str, RoleSpec()]	The path of the lattice of oscillators.	required

Examples:

>>> from skcomplex.systems.biology.neuroscience import InhibitorDriveInteraction
>>> InhibitorDriveInteraction('drive', inhibitor='z', oscillators='grid').reduce
'max'

Source code in src/skcomplex/systems/biology/neuroscience/_networks.py

def __init__(
    self,
    label: str,
    *,
    inhibitor: Annotated[str, RoleSpec()],
    oscillators: Annotated[str, RoleSpec()],
) -> None:
    super().__init__(
        label,
        target=inhibitor,
        source=oscillators,
        write='z',
        function=self._get_drive,
        reduce='max',
        mode='group',
    )
    self.inhibitor = inhibitor
    self.oscillators = oscillators

InhibitorElement(label, visual=None, *, z=0.0, phi=3.0)

Bases: Element

Global inhibitor of a LEGION network, whose activity z decays at rate phi without drive.

Parameters:

Name	Type	Description	Default
label	str	The label of the inhibitor.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
z	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial activity.	0.0
phi	Annotated[float, ParameterSpec(domain=open(0, oo))]	The rate of the inhibitor.	3.0

Examples:

>>> from skcomplex.systems.biology.neuroscience import InhibitorElement
>>> InhibitorElement('inhibitor').phi
3.0

Source code in src/skcomplex/systems/biology/neuroscience/_networks.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    z: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    phi: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 3.0,
) -> None:
    super().__init__(label, visual=visual)
    self.z = z
    self.phi = phi

IzhikevichElement(label, visual=None, *, v=-65.0, u=-13.0, a=0.02, b=0.2, I=10.0)

Bases: Element

Izhikevich neuron, a quadratic potential v and a linear recovery u.

The potential follows v' = 0.04 v² + 5 v + 140 - u + I and the recovery u' = a (b v - u). A spike reset fires the neuron.

Parameters:

Name	Type	Description	Default
label	str	The label of the neuron.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
v	Annotated[float, StateSpec()]	The initial potential.	-65.0
u	Annotated[float, StateSpec()]	The initial recovery.	-13.0
a	Annotated[float, ParameterSpec()]	The time scale of the recovery.	0.02
b	Annotated[float, ParameterSpec()]	The sensitivity of the recovery to the potential.	0.2
I	Annotated[float, ParameterSpec()]	The input current.	10.0

Examples:

>>> from skcomplex.systems.biology.neuroscience import IzhikevichElement
>>> IzhikevichElement('neuron').b
0.2

Source code in src/skcomplex/systems/biology/neuroscience/_neurons.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    v: Annotated[float, StateSpec()] = -65.0,
    u: Annotated[float, StateSpec()] = -13.0,
    a: Annotated[float, ParameterSpec()] = 0.02,
    b: Annotated[float, ParameterSpec()] = 0.2,
    I: Annotated[float, ParameterSpec()] = 10.0,
) -> None:
    super().__init__(label, visual=visual)
    self.v = v
    self.u = u
    self.a = a
    self.b = b
    self.I = I

LeakyIntegrateAndFireElement(label, visual=None, *, v=-65.0, rest=-65.0, resistance=10.0, tau=10.0, I=2.0)

Bases: Element

Leaky integrate-and-fire neuron, whose potential relaxes to rest + resistance I over the time tau.

The potential follows tau v' = -(v - rest) + resistance I. A spike reset fires the neuron. Units are mV, ms, MΩ, and nA.

Parameters:

Name	Type	Description	Default
label	str	The label of the neuron.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
v	Annotated[float, StateSpec()]	The initial potential.	-65.0
rest	Annotated[float, ParameterSpec()]	The resting potential.	-65.0
resistance	Annotated[float, ParameterSpec(domain=open(0, oo))]	The membrane resistance.	10.0
tau	Annotated[float, ParameterSpec(domain=open(0, oo))]	The membrane time constant.	10.0
I	Annotated[float, ParameterSpec()]	The input current.	2.0

Examples:

>>> from skcomplex.systems.biology.neuroscience import LeakyIntegrateAndFireElement
>>> LeakyIntegrateAndFireElement('neuron').tau
10.0

Source code in src/skcomplex/systems/biology/neuroscience/_neurons.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    v: Annotated[float, StateSpec()] = -65.0,
    rest: Annotated[float, ParameterSpec()] = -65.0,
    resistance: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 10.0,
    tau: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 10.0,
    I: Annotated[float, ParameterSpec()] = 2.0,
) -> None:
    super().__init__(label, visual=visual)
    self.v = v
    self.rest = rest
    self.resistance = resistance
    self.tau = tau
    self.I = I

MutualInhibitionInteraction(label, *, units, beta=2.0)

Bases: CouplingInteraction

Inhibition of every unit by the summed activity of all the other units, -beta x_i sum_{j != i} x_j / tau_i.

It is the coupling of the units to each other, as separate elements or as the members of a group, with the term -beta x_i x_j / tau_i of every other unit j.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
units	Annotated[str | Sequence[str], RoleSpec()]	The paths of the units, an element set, or the path of a group of units.	required
beta	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The strength of the inhibition.	2.0

Examples:

>>> from skcomplex.systems.biology.neuroscience import MutualInhibitionInteraction
>>> MutualInhibitionInteraction('inhibition', units=['a', 'b'], beta=2.0).beta
2.0

Source code in src/skcomplex/systems/biology/neuroscience/_populations.py

def __init__(
    self,
    label: str,
    *,
    units: Annotated[str | Sequence[str], RoleSpec()],
    beta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 2.0,
) -> None:
    super().__init__(label, target=units, function=self._get_inhibition, mode='group')
    self.units = units
    self.beta = beta

NeuralFieldSiteElement(label, visual=None, *, u=0.0, tau=1.0, h=-0.2, A=1.0, sigma=1.0)

Bases: Element

Site of a neural field, the mean membrane potential of the neurons at one position.

The site holds the potential u, its time constant tau, the resting level h, and the amplitude A and the range sigma of the lateral-inhibition kernel w(x) = A (1 - |x| / sigma) exp(-|x| / sigma), which excites near neighbours and inhibits distant ones.

Parameters:

Name	Type	Description	Default
label	str	The label of the site or the field.	required
visual	VisualSpec | None	How the site is drawn, or None for default_visual_.	None
u	Annotated[float, StateSpec()]	The initial membrane potential.	0.0
tau	Annotated[float, ParameterSpec(domain=open(0, oo))]	The time constant of the potential.	1.0
h	Annotated[float, ParameterSpec()]	The resting level.	-0.2
A	Annotated[float, ParameterSpec(domain=open(0, oo))]	The amplitude of the kernel.	1.0
sigma	Annotated[float, ParameterSpec(domain=open(0, oo))]	The range of the kernel.	1.0

Examples:

>>> from skcomplex.systems.biology.neuroscience import NeuralFieldSiteElement
>>> NeuralFieldSiteElement('site').h
-0.2

Source code in src/skcomplex/systems/biology/neuroscience/_neural_fields.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    u: Annotated[float, StateSpec()] = 0.0,
    tau: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    h: Annotated[float, ParameterSpec()] = -0.2,
    A: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    sigma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.u = u
    self.tau = tau
    self.h = h
    self.A = A
    self.sigma = sigma

NeuralMassElement(label, visual=None, *, r=0.01, v=-1.0, s=0.0, tau=20.0, delta=1.0, H=-5.0, tau_d=10.0)

Bases: Element

Population of quadratic integrate-and-fire neurons in the exact mean field of Montbrió, Pazó, and Roxin.

The excitabilities of the neurons follow a Lorentzian distribution of median H and half-width delta, so the firing rate r and the mean membrane potential v of the population follow dr/dt = delta / (tau^2 pi) + 2 r v / tau and dv/dt = (v^2 + H) / tau - tau (pi r)^2 exactly, with the membrane time constant tau. The synaptic field s filters the rate with exponential synapses, ds/dt = (r - s) / tau_d. Couplings and inputs add to the rate of v.

Parameters:

Name	Type	Description	Default
label	str	The label of the population.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
r	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial firing rate, in spikes per millisecond.	0.01
v	Annotated[float, StateSpec()]	The initial mean membrane potential.	-1.0
s	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial synaptic field.	0.0
tau	Annotated[float, ParameterSpec(domain=open(0, oo))]	The membrane time constant, in milliseconds.	20.0
delta	Annotated[float, ParameterSpec(domain=open(0, oo))]	The half-width of the Lorentzian excitabilities.	1.0
H	Annotated[float, ParameterSpec()]	The median of the Lorentzian excitabilities.	-5.0
tau_d	Annotated[float, ParameterSpec(domain=open(0, oo))]	The decay time of the synapses, in milliseconds.	10.0

Examples:

>>> from skcomplex.systems.biology.neuroscience import NeuralMassElement
>>> NeuralMassElement('excitatory', tau=20.0).tau
20.0

Source code in src/skcomplex/systems/biology/neuroscience/_neural_masses.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    r: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.01,
    v: Annotated[float, StateSpec()] = -1.0,
    s: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
    tau: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 20.0,
    delta: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    H: Annotated[float, ParameterSpec()] = -5.0,
    tau_d: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 10.0,
) -> None:
    super().__init__(label, visual=visual)
    self.r = r
    self.v = v
    self.s = s
    self.tau = tau
    self.delta = delta
    self.H = H
    self.tau_d = tau_d

PhaseUnitElement(label, visual=None, *, theta=0.0, omega=1.0)

Bases: Element

Neural unit reduced to a phase that turns at its natural frequency, theta' = omega.

Couplings to other units add to the rate of the phase.

Parameters:

Name	Type	Description	Default
label	str	The label of the unit or of the group of units.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
theta	Annotated[float, StateSpec(period=2 * pi)]	The initial phase.	0.0
omega	Annotated[float, ParameterSpec()]	The natural frequency.	1.0

Examples:

>>> from skcomplex.systems.biology.neuroscience import PhaseUnitElement
>>> PhaseUnitElement('unit', omega=1.5).omega
1.5

Source code in src/skcomplex/systems/biology/neuroscience/_plastic_networks.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    theta: Annotated[float, StateSpec(period=2 * pi)] = 0.0,
    omega: Annotated[float, ParameterSpec()] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.theta = theta
    self.omega = omega

PolyhomeostaticNeuronElement(label, visual=None, *, x=0.0, a=1.0, b=0.0, mu=0.2, eps_a=0.05, eps_b=0.05, Gamma=1.0, I0=0.0, A=1.0, omega=1.0)

Bases: Element

Rate neuron that adapts its gain and threshold so that its firing rate becomes exponentially distributed.

The membrane potential x integrates the input I0 + A sin(omega t) with rate Gamma, and the neuron fires at rate y = 1 / (1 + exp(a (b - x))). Stochastic gradient descent on the Kullback-Leibler divergence between the distribution of y and an exponential of mean mu gives a' = eps_a (1/a + (x - b) d) and b' = -eps_b a d with d = 1 - (2 + 1/mu) y + y²/mu. The element holds the states and parameters and has no laws of its own: the relaxation and plasticity interactions of polyhomeostatic_neuron read its parameters.

Parameters:

Name	Type	Description	Default
label	str	The label of the neuron.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec()]	The initial membrane potential.	0.0
a	Annotated[float, StateSpec(domain=open(0, oo))]	The initial gain.	1.0
b	Annotated[float, StateSpec()]	The initial threshold.	0.0
mu	Annotated[float, ParameterSpec(domain=open(0, oo))]	The target mean firing rate.	0.2
eps_a	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The adaptation rate of the gain.	0.05
eps_b	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The adaptation rate of the threshold.	0.05
Gamma	Annotated[float, ParameterSpec(domain=open(0, oo))]	The leak rate of the membrane potential.	1.0
I0	Annotated[float, ParameterSpec()]	The mean input.	0.0
A	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The amplitude of the periodic input.	1.0
omega	Annotated[float, ParameterSpec(domain=open(0, oo))]	The angular frequency of the periodic input.	1.0

Examples:

>>> from skcomplex.systems.biology.neuroscience import PolyhomeostaticNeuronElement
>>> PolyhomeostaticNeuronElement('neuron').mu
0.2

Source code in src/skcomplex/systems/biology/neuroscience/_populations.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec()] = 0.0,
    a: Annotated[float, StateSpec(domain=RealSet.open(0, oo))] = 1.0,
    b: Annotated[float, StateSpec()] = 0.0,
    mu: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.2,
    eps_a: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.05,
    eps_b: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.05,
    Gamma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    I0: Annotated[float, ParameterSpec()] = 0.0,
    A: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    omega: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.a = a
    self.b = b
    self.mu = mu
    self.eps_a = eps_a
    self.eps_b = eps_b
    self.Gamma = Gamma
    self.I0 = I0
    self.A = A
    self.omega = omega

PulseCouplingInteraction(label, *, neurons, kappa=1.0, n=2)

Bases: Interaction

Global coupling of theta neurons through the mean of the smooth pulses that every neuron emits.

Every neuron receives the current (1 + cos theta) kappa I, where I = mean_j a_n (1 - cos theta_j)^n is the mean pulse of all neurons at the current state. The pulse peaks when a neuron fires, and a_n gives it a unit mean over the phase. With kappa > 0 the coupling is excitatory.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
neurons	Annotated[str, RoleSpec()]	The path of the group of theta neurons.	required
kappa	Annotated[float, ParameterSpec()]	The coupling strength.	1.0
n	int	The sharpness of the pulses, a positive integer.	2

Examples:

>>> from skcomplex.systems.biology.neuroscience import PulseCouplingInteraction
>>> PulseCouplingInteraction('coupling', neurons='neurons', kappa=2.0).kappa
2.0

Source code in src/skcomplex/systems/biology/neuroscience/_populations.py

def __init__(
    self,
    label: str,
    *,
    neurons: Annotated[str, RoleSpec()],
    kappa: Annotated[float, ParameterSpec()] = 1.0,
    n: int = 2,
) -> None:
    super().__init__(label)
    self.neurons = neurons
    self.kappa = kappa
    self.n = n

RateUnitElement(label, visual=None, *, x=0.1, I=1.0, tau=1.0)

Bases: Element

Neural unit of a competitive network, whose activity grows with its input and saturates, tau x' = x (I - x).

Parameters:

Name	Type	Description	Default
label	str	The label of the unit.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial activity.	0.1
I	Annotated[float, ParameterSpec()]	The input.	1.0
tau	Annotated[float, ParameterSpec(domain=open(0, oo))]	The time constant.	1.0

Examples:

>>> from skcomplex.systems.biology.neuroscience import RateUnitElement
>>> RateUnitElement('unit', I=1.2).I
1.2

Source code in src/skcomplex/systems/biology/neuroscience/_populations.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
    I: Annotated[float, ParameterSpec()] = 1.0,
    tau: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.I = I
    self.tau = tau

SpikeResetInteraction(label, *, target, state='v', threshold=1.0, after=0.0, kicked=None, kick=0.0)

Bases: 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.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction and of its events.	required
target	Annotated[str, RoleSpec()]	The path of the neuron.	required
state	str	The name of the potential.	'v'
threshold	Annotated[float, ParameterSpec()]	The potential of a spike.	1.0
after	Annotated[float, ParameterSpec()]	The potential after a spike.	0.0
kicked	str | None	The name of the state a spike increments, or None for none.	None
kick	Annotated[float, ParameterSpec()]	The increment of the kicked state at a spike.	0.0

Raises:

Type	Description
ValidationError	The direction is not -1, 0, or 1, or the guard or a new value is not a function.

Examples:

>>> from skcomplex.core.components.elements import CompartmentElement
>>> from skcomplex.core.components.interactions import RateInteraction
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.runtime import Simulation
>>> from skcomplex.systems.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).prepare().run(until=2.1).trajectory_
>>> [round(event.time, 6) for event in trajectory.events]
[1.0, 2.0]

Source code in src/skcomplex/systems/biology/common/excitable/_spike_reset.py

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

SpikeTimingPlasticityInteraction(label, *, neurons, synapses, k=1.0, epsilon=0.01, A_plus=1.0, A_minus=1.05, kappa_max=1.0)

Bases: 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.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
neurons	Annotated[str, RoleSpec()]	The path of the group of traced theta neurons on the nodes.	required
synapses	Annotated[str, RoleSpec()]	The path of the group of synapses on the edges.	required
k	Annotated[float, ParameterSpec(domain=open(0, oo))]	The mean in-degree that normalizes the coupling.	1.0
epsilon	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The learning rate.	0.01
A_plus	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The amplitude of potentiation.	1.0
A_minus	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The amplitude of depression.	1.05
kappa_max	Annotated[float, ParameterSpec(domain=open(0, oo))]	The largest efficacy.	1.0

Examples:

>>> from skcomplex.systems.biology.neuroscience import SpikeTimingPlasticityInteraction
>>> SpikeTimingPlasticityInteraction('stdp', neurons='neurons', synapses='synapses').A_minus
1.05

Source code in src/skcomplex/systems/biology/neuroscience/_plastic_networks.py

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

SpikingNeuronElement(label, visual=None, *, v=0.0, s=0.0, tau=20.0, mu=25.0, tau_s=2.0)

Bases: 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.

Parameters:

Name	Type	Description	Default
label	str	The label of the neuron.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
v	Annotated[float, StateSpec()]	The initial potential.	0.0
s	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial synaptic trace.	0.0
tau	Annotated[float, ParameterSpec(domain=open(0, oo))]	The membrane time constant.	20.0
mu	Annotated[float, ParameterSpec()]	The constant external drive, the potential the membrane relaxes to.	25.0
tau_s	Annotated[float, ParameterSpec(domain=open(0, oo))]	The decay time of the synaptic trace.	2.0

Examples:

>>> from skcomplex.systems.biology.neuroscience import SpikingNeuronElement
>>> SpikingNeuronElement('neuron').tau_s
2.0

Source code in src/skcomplex/systems/biology/neuroscience/_networks.py

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

SynapseElement(label, visual=None, *, kappa=0.5)

Bases: Element

Synapse on an edge of a network, with an efficacy that plasticity changes.

Parameters:

Name	Type	Description	Default
label	str	The label of the synapse or of the group on the edges.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
kappa	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial efficacy.	0.5

Examples:

>>> from skcomplex.systems.biology.neuroscience import SynapseElement
>>> SynapseElement('synapse', kappa=0.5).kappa
0.5

Source code in src/skcomplex/systems/biology/neuroscience/_plastic_networks.py

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

SynapticDriveInteraction(label, *, target, source, J=1.0, field='r')

Bases: Interaction

Synaptic drive of a population by the rate or the synaptic field of another, dv/dt += J x.

Instantaneous synapses pass the firing rate r of the source, and exponential synapses its synaptic field s. A negative strength J is inhibition.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
target	Annotated[str, RoleSpec()]	The path of the driven population.	required
source	Annotated[str, RoleSpec()]	The path of the driving population, which may be the target.	required
J	Annotated[float, ParameterSpec()]	The strength of the synapses, negative for inhibition.	1.0
field	str	'r' for instantaneous synapses or 's' for exponential synapses.	'r'

Raises:

Type	Description
ValidationError	The field is unknown.

Examples:

>>> from skcomplex.systems.biology.neuroscience import SynapticDriveInteraction
>>> SynapticDriveInteraction('ie', target='e', source='i', J=-10.0).J
-10.0

Source code in src/skcomplex/systems/biology/neuroscience/_neural_masses.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    source: Annotated[str, RoleSpec()],
    J: Annotated[float, ParameterSpec()] = 1.0,
    field: str = 'r',
) -> None:
    super().__init__(label)
    self.target = target
    self.source = source
    self.J = J
    self.field = field

TermanWangElement(label, visual=None, *, x=SILENT, y=0.0, I=0.2, m=1.0, W=6.0, epsilon=0.02, gamma=6.0, beta=0.1, theta_x=-0.5, kappa=0.02)

Bases: Element

Relaxation oscillator of Terman and Wang, the dynamics of one site of a LEGION lattice.

The fast variable follows x' = 3x - x³ + 2 - y + I and the slow one y' = epsilon (gamma (1 + tanh(x / beta)) - y). A positive input I makes the oscillator oscillate and a negative one makes it excitable. The parameters m, W, theta_x, and kappa belong to the coupling of a lattice of oscillators, a CouplingInteraction that adds W m (sum_k m_k H(x_k - theta_x)) over the four nearest sites, where H is a smooth step of width kappa and m marks the stimulated sites, so only stimulated neighbours excite each other.

Parameters:

Name	Type	Description	Default
label	str	The label of the oscillator or the lattice field.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec()]	The initial fast variable.	SILENT
y	Annotated[float, StateSpec()]	The initial slow variable.	0.0
I	Annotated[float, ParameterSpec()]	The external input.	0.2
m	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	1 for a stimulated site and 0 otherwise.	1.0
W	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The weight of each excitatory link into the site.	6.0
epsilon	Annotated[float, ParameterSpec(domain=open(0, oo))]	The time scale of the slow variable.	0.02
gamma	Annotated[float, ParameterSpec(domain=open(0, oo))]	The height of the slow nullcline.	6.0
beta	Annotated[float, ParameterSpec(domain=open(0, oo))]	The width of the slow nullcline's step.	0.1
theta_x	Annotated[float, ParameterSpec()]	The threshold of the active phase.	-0.5
kappa	Annotated[float, ParameterSpec(domain=open(0, oo))]	The width of the smooth steps.	0.02

Examples:

>>> from skcomplex.systems.biology.neuroscience import TermanWangElement
>>> TermanWangElement('oscillator').gamma
6.0

Source code in src/skcomplex/systems/biology/neuroscience/_networks.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec()] = SILENT,
    y: Annotated[float, StateSpec()] = 0.0,
    I: Annotated[float, ParameterSpec()] = 0.2,
    m: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    W: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 6.0,
    epsilon: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.02,
    gamma: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 6.0,
    beta: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.1,
    theta_x: Annotated[float, ParameterSpec()] = -0.5,
    kappa: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.02,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.y = y
    self.I = I
    self.m = m
    self.W = W
    self.epsilon = epsilon
    self.gamma = gamma
    self.beta = beta
    self.theta_x = theta_x
    self.kappa = kappa

ThetaForcingInteraction(label, *, target, I0=10.0, nu=5.0)

Bases: Interaction

Theta forcing of a population by the current I0 (1 - cos(2 pi nu t)) / 2, as in optogenetic stimulation.

The current enters the potential of the population as its excitabilities do, dv/dt += I(t) / tau. Time is in milliseconds and the frequency nu in hertz.

Parameters:

Name	Type	Description	Default
label	str	The label of the interaction.	required
target	Annotated[str, RoleSpec()]	The path of the forced population.	required
I0	Annotated[float, ParameterSpec()]	The peak current.	10.0
nu	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The frequency of the forcing, in hertz.	5.0

Examples:

>>> from skcomplex.systems.biology.neuroscience import ThetaForcingInteraction
>>> ThetaForcingInteraction('theta', target='e', nu=8.0).nu
8.0

Source code in src/skcomplex/systems/biology/neuroscience/_neural_masses.py

def __init__(
    self,
    label: str,
    *,
    target: Annotated[str, RoleSpec()],
    I0: Annotated[float, ParameterSpec()] = 10.0,
    nu: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 5.0,
) -> None:
    super().__init__(label)
    self.target = target
    self.I0 = I0
    self.nu = nu

ThetaNeuronElement(label, visual=None, *, theta=0.0, eta=1.0)

Bases: Element

Theta neuron of Ermentrout and Kopell, theta' = (1 - cos theta) + (1 + cos theta) eta, a spiking phase.

The neuron fires when its phase passes π. With eta > 0 it fires periodically with period π / sqrt(eta), and with eta < 0 it is excitable. The change of variable V = tan(theta / 2) turns it into the quadratic integrate-and-fire neuron V' = V² + eta.

Parameters:

Name	Type	Description	Default
label	str	The label of the neuron or the group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
theta	Annotated[float, StateSpec(period=2 * pi)]	The initial phase.	0.0
eta	Annotated[float, ParameterSpec()]	The excitability.	1.0

Examples:

>>> from skcomplex.core.components.elements import GroupElement
>>> from skcomplex.systems.biology.neuroscience import ThetaNeuronElement
>>> GroupElement('neurons', element=ThetaNeuronElement, size=10, values={'eta': 0.5}).size
10

Source code in src/skcomplex/systems/biology/common/excitable/_theta_neuron.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    theta: Annotated[float, StateSpec(period=2 * pi)] = 0.0,
    eta: Annotated[float, ParameterSpec()] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.theta = theta
    self.eta = eta

TracedThetaNeuronElement(label, visual=None, *, theta=0.0, a=0.0, b=0.0, eta=1.0, tau_plus=2.0, tau_minus=2.0)

Bases: 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.

Parameters:

Name	Type	Description	Default
label	str	The label of the neuron or the group.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
theta	Annotated[float, StateSpec(period=2 * pi)]	The initial phase.	0.0
a	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial presynaptic trace.	0.0
b	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial postsynaptic trace.	0.0
eta	Annotated[float, ParameterSpec()]	The excitability.	1.0
tau_plus	Annotated[float, ParameterSpec(domain=open(0, oo))]	The time of the potentiation window.	2.0
tau_minus	Annotated[float, ParameterSpec(domain=open(0, oo))]	The time of the depression window.	2.0

Examples:

>>> from skcomplex.systems.biology.neuroscience import TracedThetaNeuronElement
>>> TracedThetaNeuronElement('neuron').tau_plus
2.0

Source code in src/skcomplex/systems/biology/neuroscience/_plastic_networks.py

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

WilsonCowanElement(label, visual=None, *, E=0.1, I=0.05, c1=16.0, c2=12.0, c3=15.0, c4=3.0, a_e=1.3, theta_e=4.0, a_i=2.0, theta_i=3.7, r_e=1.0, r_i=1.0, P=1.25, Q=0.0, tau_e=1.0, tau_i=1.0)

Bases: Element

Excitatory and inhibitory populations of Wilson and Cowan, as fractions of active cells.

The excitatory activity follows tau_e E' = -E + (k_e - r_e E) S_e(c1 E - c2 I + P) and the inhibitory activity tau_i I' = -I + (k_i - r_i I) S_i(c3 E - c4 I + Q), where S(x) = 1 / (1 + exp(-a (x - theta))) - 1 / (1 + exp(a theta)) is a sigmoid shifted to vanish at zero input and k is its largest value. The element holds the states and parameters and has no laws of its own: the relaxations of wilson_cowan read its parameters.

Parameters:

Name	Type	Description	Default
label	str	The label of the populations.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
E	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial excitatory activity.	0.1
I	Annotated[float, StateSpec(domain=closed_open(0, oo))]	The initial inhibitory activity.	0.05
c1	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The excitatory to excitatory coupling.	16.0
c2	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The inhibitory to excitatory coupling.	12.0
c3	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The excitatory to inhibitory coupling.	15.0
c4	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The inhibitory to inhibitory coupling.	3.0
a_e	Annotated[float, ParameterSpec(domain=open(0, oo))]	The gain of the excitatory response.	1.3
theta_e	Annotated[float, ParameterSpec()]	The threshold of the excitatory response.	4.0
a_i	Annotated[float, ParameterSpec(domain=open(0, oo))]	The gain of the inhibitory response.	2.0
theta_i	Annotated[float, ParameterSpec()]	The threshold of the inhibitory response.	3.7
r_e	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The refractory period of excitatory cells.	1.0
r_i	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The refractory period of inhibitory cells.	1.0
P	Annotated[float, ParameterSpec()]	The external input to the excitatory population.	1.25
Q	Annotated[float, ParameterSpec()]	The external input to the inhibitory population.	0.0
tau_e	Annotated[float, ParameterSpec(domain=open(0, oo))]	The time constant of the excitatory population.	1.0
tau_i	Annotated[float, ParameterSpec(domain=open(0, oo))]	The time constant of the inhibitory population.	1.0

Examples:

>>> from skcomplex.systems.biology.neuroscience import WilsonCowanElement
>>> WilsonCowanElement('cortex', P=1.25).c1
16.0

Source code in src/skcomplex/systems/biology/neuroscience/_populations.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    E: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.1,
    I: Annotated[float, StateSpec(domain=RealSet.closed_open(0, oo))] = 0.05,
    c1: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 16.0,
    c2: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 12.0,
    c3: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 15.0,
    c4: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 3.0,
    a_e: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.3,
    theta_e: Annotated[float, ParameterSpec()] = 4.0,
    a_i: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 2.0,
    theta_i: Annotated[float, ParameterSpec()] = 3.7,
    r_e: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    r_i: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 1.0,
    P: Annotated[float, ParameterSpec()] = 1.25,
    Q: Annotated[float, ParameterSpec()] = 0.0,
    tau_e: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    tau_i: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.E = E
    self.I = I
    self.c1 = c1
    self.c2 = c2
    self.c3 = c3
    self.c4 = c4
    self.a_e = a_e
    self.theta_e = theta_e
    self.a_i = a_i
    self.theta_i = theta_i
    self.r_e = r_e
    self.r_i = r_i
    self.P = P
    self.Q = Q
    self.tau_e = tau_e
    self.tau_i = tau_i

adaptive_delay_network(n_units=5, *, omega=1.0, spread=0.1, kappa=1.0, epsilon=0.1, rho=1.0, tau0=3.0, seed=0)

Return phase units coupled all to all through conduction delays that adapt to their activity, as myelination.

Unit i follows dtheta_i/dt = omega_i + (kappa / N) sum_j sin(theta_j(t - tau_ij) - theta_i), and every delay dtau_ij/dt = epsilon H(tau_ij) (-(tau_ij - tau0) + rho sin(theta_j - theta_i)) relaxes to its baseline tau0 and is pushed by the phase difference of its units. With the defaults, fixed delays of 3 leave the units incoherent, and the adaptation of the delays restores a partial synchrony.

Parameters:

Name	Type	Description	Default
n_units	int	The number of units.	5
omega	float	The mean natural frequency.	1.0
spread	float	The standard deviation of the natural frequencies.	0.1
kappa	float	The coupling strength, divided by the number of units.	1.0
epsilon	float	The adaptation rate of the delays.	0.1
rho	float	The strength of the push of the phase difference on the delays.	1.0
tau0	float	The baseline delay of every edge.	3.0
seed	int | None	The seed of the frequencies and the initial phases.	0

Returns:

Name	Type	Description
system	System	The system, with the phases unit_<i>.theta of the units and the delays axon_<i>_<j>.tau of the axons from unit i to unit j.

Examples:

>>> from skcomplex.systems.biology.neuroscience import adaptive_delay_network
>>> len(adaptive_delay_network(4).interactions)
12

Source code in src/skcomplex/systems/biology/neuroscience/_plastic_networks.py

@register_system(
    field='biology/neuroscience',
    tags=('adaptive delays', 'myelination', 'phase oscillators', 'adaptive network', 'dde'),
    references=('Park and Lefebvre 2020', 'Berner et al. 2023'),
    defaults={'until': 300.0, 'dt': 0.5},
)
def adaptive_delay_network(
    n_units: int = 5,
    *,
    omega: float = 1.0,
    spread: float = 0.1,
    kappa: float = 1.0,
    epsilon: float = 0.1,
    rho: float = 1.0,
    tau0: float = 3.0,
    seed: int | None = 0,
) -> System:
    """Return phase units coupled all to all through conduction delays that adapt to their activity, as myelination.

    Unit `i` follows `dtheta_i/dt = omega_i + (kappa / N) sum_j sin(theta_j(t - tau_ij) - theta_i)`, and every delay
    `dtau_ij/dt = epsilon H(tau_ij) (-(tau_ij - tau0) + rho sin(theta_j - theta_i))` relaxes to its baseline `tau0`
    and is pushed by the phase difference of its units. With the defaults, fixed delays of 3 leave the units
    incoherent, and the adaptation of the delays restores a partial synchrony.

    Args:
        n_units:
            The number of units.

        omega:
            The mean natural frequency.

        spread:
            The standard deviation of the natural frequencies.

        kappa:
            The coupling strength, divided by the number of units.

        epsilon:
            The adaptation rate of the delays.

        rho:
            The strength of the push of the phase difference on the delays.

        tau0:
            The baseline delay of every edge.

        seed:
            The seed of the frequencies and the initial phases.

    Returns:
        system:
            The system, with the phases `unit_<i>.theta` of the units and the delays `axon_<i>_<j>.tau` of the axons
            from unit `i` to unit `j`.

    Examples:
        >>> from skcomplex.systems.biology.neuroscience import adaptive_delay_network
        >>> len(adaptive_delay_network(4).interactions)
        12
    """
    rng = np.random.default_rng(seed)
    phases, frequencies = rng.uniform(0, 2 * np.pi, n_units), rng.normal(omega, spread, n_units)
    units = [
        PhaseUnitElement(f'unit_{index}', theta=float(phases[index]), omega=float(frequencies[index]))
        for index in range(n_units)
    ]
    pairs = list(permutations(range(n_units), 2))
    axons = [AxonElement(f'axon_{source}_{target}', tau=tau0, tau0=tau0) for source, target in pairs]
    adaptations = [
        DelayAdaptationInteraction(
            f'myelination_{source}_{target}',
            source=f'unit_{source}',
            target=f'unit_{target}',
            axon=f'axon_{source}_{target}',
            kappa=kappa / n_units,
            epsilon=epsilon,
            rho=rho,
        )
        for source, target in pairs
    ]
    return System(elements=[*units, *axons], interactions=adaptations)

adaptrode(*, alpha0=0.25, alpha1=0.05, alpha2=0.01, delta0=0.1, delta1=0.02, delta2=0.002, period=10.0, stop=150.0)

Return an Adaptrode synapse driven by a train of input pulses that stops after a while.

While the pulses arrive the fast weight follows each pulse, the medium weight builds up, and the slow weight builds up further. After the input stops the fast weight falls back quickly to the medium one, which keeps a trace of the activity far longer, while the slow weight retains a long-term memory.

Parameters:

Name	Type	Description	Default
alpha0	float	The learning rate of the fast weight.	0.25
alpha1	float	The learning rate of the medium weight.	0.05
alpha2	float	The learning rate of the slow weight.	0.01
delta0	float	The decay rate of the fast weight.	0.1
delta1	float	The decay rate of the medium weight.	0.02
delta2	float	The decay rate of the slow weight.	0.002
period	float	The number of steps between the starts of two pulses.	10.0
stop	float	The step after which the input stays off.	150.0

Returns:

Name	Type	Description
system	System	The system, with the weights synapse.w0, synapse.w1, and synapse.w2.

Examples:

>>> from skcomplex.systems.biology.neuroscience import adaptrode
>>> adaptrode(stop=50.0).elements[0].stop
50.0

Source code in src/skcomplex/systems/biology/neuroscience/_synapses.py

@register_system(
    field='biology/neuroscience',
    tags=('synapse', 'learning', 'memory', 'discrete map'),
    references=('Mobus 1994', 'Mobus 2022'),
    defaults={'until': 400, 'dt': 1},
)
def adaptrode(
    *,
    alpha0: float = 0.25,
    alpha1: float = 0.05,
    alpha2: float = 0.01,
    delta0: float = 0.1,
    delta1: float = 0.02,
    delta2: float = 0.002,
    period: float = 10.0,
    stop: float = 150.0,
) -> System:
    """Return an Adaptrode synapse driven by a train of input pulses that stops after a while.

    While the pulses arrive the fast weight follows each pulse, the medium weight builds up, and the slow weight
    builds up further. After the input stops the fast weight falls back quickly to the medium one, which keeps a
    trace of the activity far longer, while the slow weight retains a long-term memory.

    Args:
        alpha0:
            The learning rate of the fast weight.

        alpha1:
            The learning rate of the medium weight.

        alpha2:
            The learning rate of the slow weight.

        delta0:
            The decay rate of the fast weight.

        delta1:
            The decay rate of the medium weight.

        delta2:
            The decay rate of the slow weight.

        period:
            The number of steps between the starts of two pulses.

        stop:
            The step after which the input stays off.

    Returns:
        system:
            The system, with the weights `synapse.w0`, `synapse.w1`, and `synapse.w2`.

    Examples:
        >>> from skcomplex.systems.biology.neuroscience import adaptrode
        >>> adaptrode(stop=50.0).elements[0].stop
        50.0
    """
    synapse = AdaptrodeElement(
        'synapse',
        alpha0=alpha0,
        alpha1=alpha1,
        alpha2=alpha2,
        delta0=delta0,
        delta1=delta1,
        delta2=delta2,
        period=period,
        stop=stop,
    )
    levels = [
        MapInteraction(
            'fast_level',
            target='synapse',
            state='w0',
            function=lambda s: s.target.w0
            + s.target.alpha0 * _input(s) * (s.target.w_max - s.target.w0)
            - s.target.delta0 * (s.target.w0 - s.target.w1),
        ),
        MapInteraction(
            'medium_level',
            target='synapse',
            state='w1',
            function=lambda s: s.target.w1
            + s.target.alpha1 * (s.target.w0 - s.target.w1)
            - s.target.delta1 * (s.target.w1 - s.target.w2),
        ),
        MapInteraction(
            'slow_level',
            target='synapse',
            state='w2',
            function=lambda s: s.target.w2
            + s.target.alpha2 * (s.target.w1 - s.target.w2)
            - s.target.delta2 * s.target.w2,
        ),
    ]
    return System(elements=[synapse], interactions=levels)

amari_neural_field(n_sites=800, *, length=20.0, tau=1.0, h=-0.2, A=1.0, sigma=1.0, width=1.0)

Return the Amari neural field, whose lateral inhibition holds a stationary bump of activity.

The potential follows tau u_t = -u + integral w(x - y) H(u(y)) dy + h on a periodic line, where H is the Heaviside firing rate, h < 0 the resting level, and w(x) = A (1 - |x| / sigma) exp(-|x| / sigma) the lateral-inhibition kernel. Amari showed that a bump of width a is stationary when W(a) + h = 0, with W(a) = A a exp(-a / sigma) the integral of w from 0 to a, and stable when w(a) < 0. A localized excitation then relaxes to the wide root, a = 2.54 for the defaults.

Parameters:

Name	Type	Description	Default
n_sites	int	The number of sites of the line.	800
length	float	The length of the line.	20.0
tau	float	The time constant of the potential.	1.0
h	float	The resting level.	-0.2
A	float	The amplitude of the kernel.	1.0
sigma	float	The range of the kernel.	1.0
width	float	The width of the initial excitation, a potential of 1 around the middle of the line and h elsewhere.	1.0

Returns:

Name	Type	Description
system	System	The system, with the potential field.u.

Examples:

>>> from skcomplex.systems.biology.neuroscience import amari_neural_field
>>> amari_neural_field(100).elements[0].values['h']
-0.2

Source code in src/skcomplex/systems/biology/neuroscience/_neural_fields.py

@register_system(
    field='biology/neuroscience',
    tags=('neural field', 'field', 'convolution', 'bump', 'working memory', 'pde'),
    references=('Amari 1977', 'arXiv:2112.12147'),
    defaults={'until': 40.0, 'dt': 0.5},
)
def amari_neural_field(
    n_sites: int = 800,
    *,
    length: float = 20.0,
    tau: float = 1.0,
    h: float = -0.2,
    A: float = 1.0,
    sigma: float = 1.0,
    width: float = 1.0,
) -> System:
    """Return the Amari neural field, whose lateral inhibition holds a stationary bump of activity.

    The potential follows `tau u_t = -u + integral w(x - y) H(u(y)) dy + h` on a periodic line, where `H` is the
    Heaviside firing rate, `h < 0` the resting level, and `w(x) = A (1 - |x| / sigma) exp(-|x| / sigma)` the
    lateral-inhibition kernel. Amari showed that a bump of width `a` is stationary when `W(a) + h = 0`, with
    `W(a) = A a exp(-a / sigma)` the integral of `w` from 0 to `a`, and stable when `w(a) < 0`. A localized excitation
    then relaxes to the wide root, `a = 2.54` for the defaults.

    Args:
        n_sites:
            The number of sites of the line.

        length:
            The length of the line.

        tau:
            The time constant of the potential.

        h:
            The resting level.

        A:
            The amplitude of the kernel.

        sigma:
            The range of the kernel.

        width:
            The width of the initial excitation, a potential of 1 around the middle of the line and `h` elsewhere.

    Returns:
        system:
            The system, with the potential `field.u`.

    Examples:
        >>> from skcomplex.systems.biology.neuroscience import amari_neural_field
        >>> amari_neural_field(100).elements[0].values['h']
        -0.2
    """
    spacing = length / n_sites
    positions = (np.arange(n_sites) + 0.5) * spacing - length / 2
    values = {'u': np.where(np.abs(positions) < width / 2, 1.0, h), 'tau': tau, 'h': h, 'A': A, 'sigma': sigma}
    field = FieldElement(
        'field',
        element=NeuralFieldSiteElement,
        space=LatticeSpace(n_sites, spacing=spacing),
        values=values,
    )
    rate = RateInteraction('rate', target='field', state='u', rate=_amari_rate)
    return System(elements=[field], interactions=[rate])

associative_thought_process(cliques=CLIQUES, *, w=0.15, z=1.0, x_c=0.85, gamma_plus=0.005, gamma_minus=0.005, start=(0, 1))

Return the dense homogeneous associative network of Gros, whose memories become active one after the other.

Every memory is a clique of sites joined by excitatory links of strength w, and every other pair of sites inhibits each other with strength z. Without reservoirs the network is a winners-take-all network whose cliques are stable attractors. The reservoirs of an active clique deplete slowly, and once they are low the clique loses its excitation and its hold on the other sites, so the activity hops to an associated clique that shares sites with it. The memories become transient attractors, and the network runs through a sequence of them, an associative thought process. The default is the network of seven sites and six cliques of Gros (2015, figure 10.12).

Parameters:

Name	Type	Description	Default
cliques	Sequence[Sequence[int]]	The memories, each a sequence of site indices.	CLIQUES
w	float	The strength of the excitatory links within a clique.	0.15
z	float	The strength of the inhibition between sites of no common clique.	1.0
x_c	float	The activity above which a site counts as active.	0.85
gamma_plus	float	The recovery rate of the reservoirs.	0.005
gamma_minus	float	The depletion rate of the reservoirs.	0.005
start	Sequence[int]	The sites that are active at the start.	(0, 1)

Returns:

Name	Type	Description
system	System	The system, with the activities sites.x and the reservoirs sites.phi.

Examples:

>>> from skcomplex.systems.biology.neuroscience import associative_thought_process
>>> associative_thought_process().elements[0].prepare().size_
7

Source code in src/skcomplex/systems/biology/neuroscience/_associative.py

@register_system(
    field='biology/neuroscience',
    tags=('associative memory', 'transient attractors', 'winners-take-all', 'cognitive systems', 'ode'),
    references=('Gros 2007', 'Gros 2015'),
    defaults={'until': 1500.0, 'dt': 2.0},
)
def associative_thought_process(
    cliques: Sequence[Sequence[int]] = CLIQUES,
    *,
    w: float = 0.15,
    z: float = 1.0,
    x_c: float = 0.85,
    gamma_plus: float = 0.005,
    gamma_minus: float = 0.005,
    start: Sequence[int] = (0, 1),
) -> System:
    """Return the dense homogeneous associative network of Gros, whose memories become active one after the other.

    Every memory is a clique of sites joined by excitatory links of strength `w`, and every other pair of sites
    inhibits each other with strength `z`. Without reservoirs the network is a winners-take-all network whose cliques
    are stable attractors. The reservoirs of an active clique deplete slowly, and once they are low the clique loses
    its excitation and its hold on the other sites, so the activity hops to an associated clique that shares sites with
    it. The memories become transient attractors, and the network runs through a sequence of them, an associative
    thought process. The default is the network of seven sites and six cliques of Gros (2015, figure 10.12).

    Args:
        cliques:
            The memories, each a sequence of site indices.

        w:
            The strength of the excitatory links within a clique.

        z:
            The strength of the inhibition between sites of no common clique.

        x_c:
            The activity above which a site counts as active.

        gamma_plus:
            The recovery rate of the reservoirs.

        gamma_minus:
            The depletion rate of the reservoirs.

        start:
            The sites that are active at the start.

    Returns:
        system:
            The system, with the activities `sites.x` and the reservoirs `sites.phi`.

    Examples:
        >>> from skcomplex.systems.biology.neuroscience import associative_thought_process
        >>> associative_thought_process().elements[0].prepare().size_
        7
    """
    size = 1 + max(max(clique) for clique in cliques)
    links = np.zeros((size, size))
    for clique in cliques:
        for first, second in combinations(clique, 2):
            links[first, second] = links[second, first] = w
    activities = np.zeros(size)
    activities[list(start)] = 1.0
    sites = GroupElement(
        'sites',
        element=AssociativeSiteElement,
        size=size,
        values={'x': activities, 'phi': 1.0, 'x_c': x_c, 'gamma_plus': gamma_plus, 'gamma_minus': gamma_minus},
    )
    competition = CliqueCompetitionInteraction('competition', sites='sites', w=tuple(map(tuple, links.tolist())), z=z)
    return System(elements=[sites], interactions=[competition])

delayed_lif_network(n_neurons=10, *, probability=0.3, J=0.5, D=1.5, mu=25.0, random_state=0)

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.

Parameters:

Name	Type	Description	Default
n_neurons	int	The number of neurons.	10
probability	float	The probability of each directed link.	0.3
J	float	The efficacy of every link in mV.	0.5
D	float	The transmission delay in ms.	1.5
mu	float	The constant drive in mV.	25.0
random_state	int | Generator | None	The seed or generator of the links and the initial potentials.	0

Returns:

Name	Type	Description
system	System	The system, with the potentials neuron_<i>.v, the traces neuron_<i>.s, and the spikes of neuron i as spike_<i> events.

Examples:

>>> from skcomplex.systems.biology.neuroscience import delayed_lif_network
>>> len(delayed_lif_network(5).elements)
5

Source code in src/skcomplex/systems/biology/neuroscience/_networks.py

@register_system(
    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 system, with the potentials `neuron_<i>.v`, the traces `neuron_<i>.s`, and the spikes of neuron `i` as
            `spike_<i>` events.

    Examples:
        >>> from skcomplex.systems.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])

fitzhugh_nagumo_ode(v=-1.2, w=-0.6, *, a=0.7, b=0.8, epsilon=0.08, I=0.5)

Return the FitzHugh-Nagumo neuron, the two-variable reduction of Hodgkin and Huxley's excitable membrane.

Without current the neuron rests on a stable fixed point and answers a large enough kick with one spike. A current I between about 0.33 and 1.42 destabilizes the rest and the neuron fires periodically.

Parameters:

Name	Type	Description	Default
v	float	The initial potential.	-1.2
w	float	The initial recovery.	-0.6
a	float	The offset of the recovery.	0.7
b	float	The damping of the recovery.	0.8
epsilon	float	The time scale of the recovery.	0.08
I	float	The applied current.	0.5

Returns:

Name	Type	Description
system	System	The system, with the potential neuron.v and the recovery neuron.w.

Examples:

>>> from skcomplex.systems.biology.neuroscience import fitzhugh_nagumo_ode
>>> fitzhugh_nagumo_ode(I=0.0).interactions[0].I
0.0

Source code in src/skcomplex/systems/biology/neuroscience/_conductances.py

@register_system(
    field='biology/neuroscience',
    tags=('neuron', 'excitability', 'relaxation oscillations'),
    references=('FitzHugh 1961', 'Nagumo, Arimoto, and Yoshizawa 1962', 'Trzaska 2023', 'Datseris and Parlitz 2022'),
    defaults={'until': 100.0, 'dt': 0.2},
)
def fitzhugh_nagumo_ode(
    v: float = -1.2,
    w: float = -0.6,
    *,
    a: float = 0.7,
    b: float = 0.8,
    epsilon: float = 0.08,
    I: float = 0.5,
) -> System:
    """Return the FitzHugh-Nagumo neuron, the two-variable reduction of Hodgkin and Huxley's excitable membrane.

    Without current the neuron rests on a stable fixed point and answers a large enough kick with one spike. A current
    `I` between about 0.33 and 1.42 destabilizes the rest and the neuron fires periodically.

    Args:
        v:
            The initial potential.

        w:
            The initial recovery.

        a:
            The offset of the recovery.

        b:
            The damping of the recovery.

        epsilon:
            The time scale of the recovery.

        I:
            The applied current.

    Returns:
        system:
            The system, with the potential `neuron.v` and the recovery `neuron.w`.

    Examples:
        >>> from skcomplex.systems.biology.neuroscience import fitzhugh_nagumo_ode
        >>> fitzhugh_nagumo_ode(I=0.0).interactions[0].I
        0.0
    """
    kinetics = FitzHughNagumoInteraction('kinetics', target='neuron', a=a, b=b, epsilon=epsilon, I=I)
    return System(elements=[FitzHughNagumoNeuronElement('neuron', v=v, w=w)], interactions=[kinetics])

hodgkin_huxley(v=-65.0, *, I=10.0, C=1.0, g_na=120.0, g_k=36.0, g_l=0.3)

Return the Hodgkin-Huxley neuron, whose sodium and potassium channels produce action potentials.

The gates start at their steady state at the initial potential, and the current is switched on at the start. A constant current above about 6.2 µA/cm² makes the neuron fire periodically, and a weaker one gives at most a few spikes before it settles.

Parameters:

Name	Type	Description	Default
v	float	The initial membrane potential in mV.	-65.0
I	float	The applied current in µA/cm².	10.0
C	float	The membrane capacitance in µF/cm².	1.0
g_na	float	The maximal sodium conductance in mS/cm².	120.0
g_k	float	The maximal potassium conductance in mS/cm².	36.0
g_l	float	The leak conductance in mS/cm².	0.3

Returns:

Name	Type	Description
system	System	The system, with the potential neuron.v, the gates neuron.n, neuron.m, and neuron.h, and its spikes as spike events.

Examples:

>>> from skcomplex.systems.biology.neuroscience import hodgkin_huxley
>>> round(hodgkin_huxley().elements[0].h, 4)
0.5961

Source code in src/skcomplex/systems/biology/neuroscience/_conductances.py

@register_system(
    field='biology/neuroscience',
    tags=('neuron', 'excitability', 'spiking', 'conductance'),
    references=('Hodgkin and Huxley 1952', 'Rinzel and Miller 1980', 'Trzaska 2023', 'arXiv:2112.12147'),
    defaults={'until': 50.0, 'dt': 0.1},
)
def hodgkin_huxley(
    v: float = -65.0,
    *,
    I: float = 10.0,
    C: float = 1.0,
    g_na: float = 120.0,
    g_k: float = 36.0,
    g_l: float = 0.3,
) -> System:
    """Return the Hodgkin-Huxley neuron, whose sodium and potassium channels produce action potentials.

    The gates start at their steady state at the initial potential, and the current is switched on at the start. A
    constant current above about 6.2 µA/cm² makes the neuron fire periodically, and a weaker one gives at most a few
    spikes before it settles.

    Args:
        v:
            The initial membrane potential in mV.

        I:
            The applied current in µA/cm².

        C:
            The membrane capacitance in µF/cm².

        g_na:
            The maximal sodium conductance in mS/cm².

        g_k:
            The maximal potassium conductance in mS/cm².

        g_l:
            The leak conductance in mS/cm².

    Returns:
        system:
            The system, with the potential `neuron.v`, the gates `neuron.n`, `neuron.m`, and `neuron.h`, and its spikes
            as `spike` events.

    Examples:
        >>> from skcomplex.systems.biology.neuroscience import hodgkin_huxley
        >>> round(hodgkin_huxley().elements[0].h, 4)
        0.5961
    """
    gates = _resting_gates(v)
    return System(
        elements=[
            HodgkinHuxleyElement(
                'neuron',
                v=v,
                I=I,
                C=C,
                g_na=g_na,
                g_k=g_k,
                g_l=g_l,
                n=gates['n'],
                m=gates['m'],
                h=gates['h'],
            ),
        ],
        interactions=[
            RateInteraction('applied_current', target='neuron', state='v', rate=lambda s: s.target.I, divisor='C'),
            RateInteraction(
                'sodium_current',
                target='neuron',
                state='v',
                rate=lambda s: -s.target.g_na * s.target.m**3 * s.target.h * (s.target.v - s.target.e_na),
                divisor='C',
            ),
            RateInteraction(
                'potassium_current',
                target='neuron',
                state='v',
                rate=lambda s: -s.target.g_k * s.target.n**4 * (s.target.v - s.target.e_k),
                divisor='C',
            ),
            RateInteraction(
                'leak_current',
                target='neuron',
                state='v',
                rate=lambda s: -s.target.g_l * (s.target.v - s.target.e_l),
                divisor='C',
            ),
            *[
                RateInteraction(f'{gate}_gating', target='neuron', state=gate, rate=_gate_kinetics(gate))
                for gate in ('n', 'm', 'h')
            ],
            CrossingInteraction('spike', target='neuron', guard=lambda s: s.target.v, direction=1),
        ],
    )

integrate_and_fire(v=-65.0, *, rest=-65.0, threshold=-50.0, reset=-65.0, resistance=10.0, tau=10.0, I=2.0)

Return the leaky integrate-and-fire neuron, tau v' = -(v - rest) + R I, which spikes and resets at a threshold.

With R I above threshold - rest it fires periodically with period tau ln(R I / (R I - threshold + rest)) when reset = rest.

Parameters:

Name	Type	Description	Default
v	float	The initial membrane potential in mV.	-65.0
rest	float	The resting potential in mV.	-65.0
threshold	float	The potential at which the neuron spikes in mV.	-50.0
reset	float	The potential after a spike in mV.	-65.0
resistance	float	The membrane resistance in MΩ.	10.0
tau	float	The membrane time constant in ms.	10.0
I	float	The input current in nA.	2.0

Returns:

Name	Type	Description
system	System	The system, with the membrane potential neuron.v and its spikes recorded as spike events.

Examples:

>>> from skcomplex.systems.biology.neuroscience import integrate_and_fire
>>> from skcomplex.core.engine.runtime import Simulation
>>> trajectory = Simulation(integrate_and_fire(), dt=0.1).prepare().run(until=50.0).trajectory_
>>> round(trajectory.events[0].time, 2)
13.86

Source code in src/skcomplex/systems/biology/neuroscience/_neurons.py

@register_system(
    field='biology/neuroscience',
    tags=('neuron', 'spiking', 'events'),
    references=('Lapicque 1907', 'Gerstner et al. 2014'),
    defaults={'until': 100.0, 'dt': 0.1},
)
def integrate_and_fire(
    v: float = -65.0,
    *,
    rest: float = -65.0,
    threshold: float = -50.0,
    reset: float = -65.0,
    resistance: float = 10.0,
    tau: float = 10.0,
    I: float = 2.0,
) -> System:
    """Return the leaky integrate-and-fire neuron, `tau v' = -(v - rest) + R I`, which spikes and resets at a threshold.

    With `R I` above `threshold - rest` it fires periodically with period `tau ln(R I / (R I - threshold + rest))`
    when `reset = rest`.

    Args:
        v:
            The initial membrane potential in mV.

        rest:
            The resting potential in mV.

        threshold:
            The potential at which the neuron spikes in mV.

        reset:
            The potential after a spike in mV.

        resistance:
            The membrane resistance in MΩ.

        tau:
            The membrane time constant in ms.

        I:
            The input current in nA.

    Returns:
        system:
            The system, with the membrane potential `neuron.v` and its spikes recorded as `spike` events.

    Examples:
        >>> from skcomplex.systems.biology.neuroscience import integrate_and_fire
        >>> from skcomplex.core.engine.runtime import Simulation
        >>> trajectory = Simulation(integrate_and_fire(), dt=0.1).prepare().run(until=50.0).trajectory_
        >>> round(trajectory.events[0].time, 2)
        13.86
    """
    neuron = LeakyIntegrateAndFireElement('neuron', v=v, rest=rest, resistance=resistance, tau=tau, I=I)
    spike = SpikeResetInteraction('spike', target='neuron', threshold=threshold, after=reset)
    return System(elements=[neuron], interactions=[spike])

izhikevich(v=-65.0, u=None, *, a=0.02, b=0.2, c=-65.0, d=8.0, I=10.0)

Return the Izhikevich neuron, whose four parameters give regular spiking, bursting, chattering, and more.

The potential follows v' = 0.04 v² + 5 v + 140 - u + I and the recovery u' = a (b v - u). When v reaches 30 mV the neuron spikes, v resets to c, and u grows by d. The defaults give regular spiking.

Parameters:

Name	Type	Description	Default
v	float	The initial membrane potential in mV.	-65.0
u	float | None	The initial recovery variable, or None for b v.	None
a	float	The time scale of the recovery.	0.02
b	float	The sensitivity of the recovery to the potential.	0.2
c	float	The potential after a spike.	-65.0
d	float	The increase of the recovery after a spike.	8.0
I	float	The input current.	10.0

Returns:

Name	Type	Description
system	System	The system, with the potential neuron.v, the recovery neuron.u, and its spikes as spike events.

Examples:

>>> from skcomplex.systems.biology.neuroscience import izhikevich
>>> izhikevich(d=2.0).interactions[0].kick
2.0

Source code in src/skcomplex/systems/biology/neuroscience/_neurons.py

@register_system(
    field='biology/neuroscience',
    tags=('neuron', 'spiking', 'events', 'bursting'),
    references=('Izhikevich 2003',),
    defaults={'until': 200.0, 'dt': 0.1},
)
def izhikevich(
    v: float = -65.0,
    u: float | None = None,
    *,
    a: float = 0.02,
    b: float = 0.2,
    c: float = -65.0,
    d: float = 8.0,
    I: float = 10.0,
) -> System:
    """Return the Izhikevich neuron, whose four parameters give regular spiking, bursting, chattering, and more.

    The potential follows `v' = 0.04 v² + 5 v + 140 - u + I` and the recovery `u' = a (b v - u)`. When `v` reaches
    30 mV the neuron spikes, `v` resets to `c`, and `u` grows by `d`. The defaults give regular spiking.

    Args:
        v:
            The initial membrane potential in mV.

        u:
            The initial recovery variable, or `None` for `b v`.

        a:
            The time scale of the recovery.

        b:
            The sensitivity of the recovery to the potential.

        c:
            The potential after a spike.

        d:
            The increase of the recovery after a spike.

        I:
            The input current.

    Returns:
        system:
            The system, with the potential `neuron.v`, the recovery `neuron.u`, and its spikes as `spike` events.

    Examples:
        >>> from skcomplex.systems.biology.neuroscience import izhikevich
        >>> izhikevich(d=2.0).interactions[0].kick
        2.0
    """
    neuron = IzhikevichElement('neuron', v=v, u=b * v if u is None else u, a=a, b=b, I=I)
    spike = SpikeResetInteraction('spike', target='neuron', threshold=PEAK, after=c, kicked='u', kick=d)
    return System(elements=[neuron], interactions=[spike])

kinouchi_copelli(graph=None, *, sigma=1.0, stimulus=0.0, n_states=5, excited=None, random_state=0)

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.

Parameters:

Name	Type	Description	Default
graph	Graph | None	The NetworkX graph, or None for an Erdős-Rényi graph of 1000 nodes with mean degree 10.	None
sigma	float	The branching ratio.	1.0
stimulus	float	The rate r of the external stimulus per step, which excites a quiescent element with probability 1 - exp(-r).	0.0
n_states	int	The number of states of the cycle, one quiescent, one excited, and n_states - 2 refractory.	5
excited	float | None	The initial fraction of excited elements, with the others quiescent, or None for states drawn uniformly over the cycle.	None
random_state	int | Generator | None	The seed or generator of the default graph and of the initial states.	0

Returns:

Name	Type	Description
system	System	The system, 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.systems.biology.neuroscience import kinouchi_copelli
>>> kinouchi_copelli(nx.random_regular_graph(4, 20, seed=0), sigma=1.0).elements[0].values['p']
0.25

Source code in src/skcomplex/systems/biology/neuroscience/_excitable.py

@register_system(
    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 system, 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.systems.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())

legion(stimulus=LEGION_OBJECTS, *, epsilon=0.1, W=6.0, W_z=1.5, random_state=0)

Return a LEGION network, a lattice of relaxation oscillators that segments a stimulus into its objects.

Stimulated oscillators excite their stimulated neighbours, so each connected object synchronizes, and a global inhibitor lets only one object be active at a time, so different objects take turns. The inhibitor is driven by the largest activity of the oscillators, so one active oscillator drives it fully. Unstimulated sites have a negative input and stay silent.

The default epsilon = 0.1 gives a cycle of about 60 time units, so with the default stimulus the second object first fires near t = 51 and the first near t = 70. The value 0.02 of Terman and Wang makes the slow phases five times longer, and then no object starts an active phase before t = 220.

Parameters:

Name	Type	Description	Default
stimulus	Sequence[Sequence[float]]	The lattice of 1 for stimulated and 0 for unstimulated sites, by rows. The default holds two objects.	LEGION_OBJECTS
epsilon	float	The time scale of the slow variables.	0.1
W	float	The total excitatory weight into a stimulated site, shared among its stimulated neighbours.	6.0
W_z	float	The strength of the global inhibition.	1.5
random_state	int | Generator | None	The seed or generator of the initial states.	0

Returns:

Name	Type	Description
system	System	The system, with the fast variables oscillators.x, the slow variables oscillators.y, and the inhibitor inhibitor.z.

Examples:

>>> from skcomplex.systems.biology.neuroscience import legion
>>> legion().elements[0].space.shape
(4, 4)

Source code in src/skcomplex/systems/biology/neuroscience/_networks.py

@register_system(
    field='biology/neuroscience',
    tags=('oscillators', 'synchronization', 'segmentation', 'lattice', 'field'),
    references=('Terman and Wang 1995', 'Wang and Terman 1995', 'Gros 2015'),
    defaults={'until': 150.0, 'dt': 0.5},
)
def legion(
    stimulus: Sequence[Sequence[float]] = LEGION_OBJECTS,
    *,
    epsilon: float = 0.1,
    W: float = 6.0,
    W_z: float = 1.5,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return a LEGION network, a lattice of relaxation oscillators that segments a stimulus into its objects.

    Stimulated oscillators excite their stimulated neighbours, so each connected object synchronizes, and a global
    inhibitor lets only one object be active at a time, so different objects take turns. The inhibitor is driven by the
    largest activity of the oscillators, so one active oscillator drives it fully. Unstimulated sites have a negative
    input and stay silent.

    The default `epsilon = 0.1` gives a cycle of about 60 time units, so with the default stimulus the second object
    first fires near `t = 51` and the first near `t = 70`. The value 0.02 of Terman and Wang makes the slow phases five
    times longer, and then no object starts an active phase before `t = 220`.

    Args:
        stimulus:
            The lattice of 1 for stimulated and 0 for unstimulated sites, by rows. The default holds two objects.

        epsilon:
            The time scale of the slow variables.

        W:
            The total excitatory weight into a stimulated site, shared among its stimulated neighbours.

        W_z:
            The strength of the global inhibition.

        random_state:
            The seed or generator of the initial states.

    Returns:
        system:
            The system, with the fast variables `oscillators.x`, the slow variables `oscillators.y`, and the inhibitor
            `inhibitor.z`.

    Examples:
        >>> from skcomplex.systems.biology.neuroscience import legion
        >>> legion().elements[0].space.shape
        (4, 4)
    """
    marks = np.asarray(stimulus, dtype=float)
    rng = np.random.default_rng(random_state)
    values = {
        'x': rng.uniform(-2.0, 2.0, marks.shape),
        'y': rng.uniform(0.0, 6.0, marks.shape),
        'I': np.where(marks > 0, 0.2, -0.02),
        'm': marks,
        'W': _link_weights(marks, W),
        'epsilon': epsilon,
    }
    oscillators = FieldElement(
        'oscillators',
        element=TermanWangElement,
        space=LatticeSpace(marks.shape),
        boundary={'x': ('fixed', SILENT), 'y': 'zero-flux'},
        values=values,
    )
    excitation = CouplingInteraction('excitation', target='oscillators', read='x', write='x', function=_excitation)
    drive = InhibitorDriveInteraction('drive', inhibitor='inhibitor', oscillators='oscillators')
    inhibition = GlobalInhibitionInteraction('inhibition', inhibitor='inhibitor', oscillators='oscillators', W_z=W_z)
    return System(
        elements=[oscillators, InhibitorElement('inhibitor')],
        interactions=[excitation, drive, inhibition],
    )

next_generation_neural_mass(setup='ping', *, H_e=1.3, H_i=None, J_ee=8.0, J_ie=10.0, J_ei=10.0, J_ii=None, delta=None, I0=None, nu=5.0)

Return the exact neural mass of quadratic integrate-and-fire neurons whose gamma rhythm nests in a theta drive.

In the 'ping' setup an excitatory and an inhibitory population, with membrane time constants of 20 and 10 ms, couple through instantaneous synapses, and the theta current drives the excitatory one. In the 'ing' setup a single inhibitory population inhibits itself through exponential synapses of 10 ms, and the current drives it. Without forcing the populations rest or oscillate collectively in the gamma band, about 45 Hz for H_e = 10, as their excitabilities cross a Hopf bifurcation. A slow theta current sweeps them across it, so gamma bursts appear at every theta peak, theta-nested gamma oscillations. The defaults are those of Segneri et al. (2020, figure 3).

Parameters:

Name	Type	Description	Default
setup	str	'ping' for excitatory and inhibitory populations or 'ing' for one inhibitory population.	'ping'
H_e	float	The median excitability of the excitatory population.	1.3
H_i	float | None	The median excitability of the inhibitory population, or None for -5 in 'ping' and 2 in 'ing'.	None
J_ee	float	The excitatory self-coupling of 'ping'.	8.0
J_ie	float	The inhibition of the excitatory population by the inhibitory one in 'ping'.	10.0
J_ei	float	The excitation of the inhibitory population by the excitatory one in 'ping'.	10.0
J_ii	float | None	The inhibitory self-coupling, or None for 0 in 'ping' and 21 in 'ing'.	None
delta	float | None	The half-width of the excitabilities, or None for 1 in 'ping' and 0.3 in 'ing'.	None
I0	float | None	The peak theta current, or None for 10 in 'ping' and 9 in 'ing'.	None
nu	float	The theta frequency, in hertz.	5.0

Returns:

Name	Type	Description
system	System	The system, with the rates excitatory.r and inhibitory.r, in spikes per millisecond, and the potentials excitatory.v and inhibitory.v, or the inhibitory population alone in 'ing'.

Raises:

Type	Description
ValidationError	The setup is unknown.

Examples:

>>> from skcomplex.systems.biology.neuroscience import next_generation_neural_mass
>>> [element.label for element in next_generation_neural_mass('ing').elements]
['inhibitory']

Source code in src/skcomplex/systems/biology/neuroscience/_neural_masses.py

@register_system(
    field='biology/neuroscience',
    tags=('neural mass', 'next generation', 'gamma oscillations', 'theta-gamma', 'ode'),
    references=('Montbrió, Pazó, and Roxin 2015', 'Segneri et al. 2020'),
    defaults={'until': 1000.0, 'dt': 0.5},
)
def next_generation_neural_mass(
    setup: str = 'ping',
    *,
    H_e: float = 1.3,
    H_i: float | None = None,
    J_ee: float = 8.0,
    J_ie: float = 10.0,
    J_ei: float = 10.0,
    J_ii: float | None = None,
    delta: float | None = None,
    I0: float | None = None,
    nu: float = 5.0,
) -> System:
    """Return the exact neural mass of quadratic integrate-and-fire neurons whose gamma rhythm nests in a theta drive.

    In the `'ping'` setup an excitatory and an inhibitory population, with membrane time constants of 20 and 10 ms,
    couple through instantaneous synapses, and the theta current drives the excitatory one. In the `'ing'` setup a
    single inhibitory population inhibits itself through exponential synapses of 10 ms, and the current drives it.
    Without forcing the populations rest or oscillate collectively in the gamma band, about 45 Hz for `H_e = 10`, as
    their excitabilities cross a Hopf bifurcation. A slow theta current sweeps them across it, so gamma bursts appear
    at every theta peak, theta-nested gamma oscillations. The defaults are those of Segneri et al. (2020, figure 3).

    Args:
        setup:
            `'ping'` for excitatory and inhibitory populations or `'ing'` for one inhibitory population.

        H_e:
            The median excitability of the excitatory population.

        H_i:
            The median excitability of the inhibitory population, or `None` for -5 in `'ping'` and 2 in `'ing'`.

        J_ee:
            The excitatory self-coupling of `'ping'`.

        J_ie:
            The inhibition of the excitatory population by the inhibitory one in `'ping'`.

        J_ei:
            The excitation of the inhibitory population by the excitatory one in `'ping'`.

        J_ii:
            The inhibitory self-coupling, or `None` for 0 in `'ping'` and 21 in `'ing'`.

        delta:
            The half-width of the excitabilities, or `None` for 1 in `'ping'` and 0.3 in `'ing'`.

        I0:
            The peak theta current, or `None` for 10 in `'ping'` and 9 in `'ing'`.

        nu:
            The theta frequency, in hertz.

    Returns:
        system:
            The system, with the rates `excitatory.r` and `inhibitory.r`, in spikes per millisecond, and the
            potentials `excitatory.v` and `inhibitory.v`, or the inhibitory population alone in `'ing'`.

    Raises:
        ValidationError:
            The setup is unknown.

    Examples:
        >>> from skcomplex.systems.biology.neuroscience import next_generation_neural_mass
        >>> [element.label for element in next_generation_neural_mass('ing').elements]
        ['inhibitory']
    """
    if setup not in SETUPS:
        error_msg = f'The setup must be one of {", ".join(SETUPS)}, not {setup!r}.'
        raise ValidationError(error_msg)
    if setup == 'ing':
        population = NeuralMassElement(
            'inhibitory',
            tau=10.0,
            delta=0.3 if delta is None else delta,
            H=2.0 if H_i is None else H_i,
        )
        return System(
            elements=[population],
            interactions=[
                SynapticDriveInteraction(
                    'ii',
                    target='inhibitory',
                    source='inhibitory',
                    J=-(21.0 if J_ii is None else J_ii),
                    field='s',
                ),
                ThetaForcingInteraction('theta', target='inhibitory', I0=9.0 if I0 is None else I0, nu=nu),
            ],
        )
    width = 1.0 if delta is None else delta
    return System(
        elements=[
            NeuralMassElement('excitatory', tau=20.0, delta=width, H=H_e),
            NeuralMassElement('inhibitory', tau=10.0, delta=width, H=-5.0 if H_i is None else H_i),
        ],
        interactions=[
            SynapticDriveInteraction('ee', target='excitatory', source='excitatory', J=J_ee),
            SynapticDriveInteraction('ie', target='excitatory', source='inhibitory', J=-J_ie),
            SynapticDriveInteraction('ei', target='inhibitory', source='excitatory', J=J_ei),
            SynapticDriveInteraction(
                'ii',
                target='inhibitory',
                source='inhibitory',
                J=-(0.0 if J_ii is None else J_ii),
            ),
            ThetaForcingInteraction('theta', target='excitatory', I0=10.0 if I0 is None else I0, nu=nu),
        ],
    )

polyhomeostatic_neuron(*, mu=0.2, eps_a=0.05, eps_b=0.05, I0=0.0, A=1.0, omega=1.0)

Return a neuron that adapts its gain and threshold so its firing rate approaches an exponential distribution.

Polyhomeostatic optimization targets the whole distribution of the firing rate, the one of maximal entropy for a given mean mu, rather than its mean alone. The neuron is driven by a periodic input, and the time average of its rate approaches mu as gain and threshold adapt.

Parameters:

Name	Type	Description	Default
mu	float	The target mean firing rate.	0.2
eps_a	float	The adaptation rate of the gain.	0.05
eps_b	float	The adaptation rate of the threshold.	0.05
I0	float	The mean input.	0.0
A	float	The amplitude of the periodic input.	1.0
omega	float	The angular frequency of the periodic input.	1.0

Returns:

Name	Type	Description
system	System	The system, with the potential neuron.x, the gain neuron.a, and the threshold neuron.b.

Examples:

>>> from skcomplex.systems.biology.neuroscience import polyhomeostatic_neuron
>>> polyhomeostatic_neuron(mu=0.1).elements[0].mu
0.1

Source code in src/skcomplex/systems/biology/neuroscience/_populations.py

@register_system(
    field='biology/neuroscience',
    tags=('neuron', 'homeostasis', 'adaptation', 'information'),
    references=('Triesch 2005', 'Markovic and Gros 2010', 'Gros 2015'),
    defaults={'until': 400.0, 'dt': 1.0},
)
def polyhomeostatic_neuron(
    *,
    mu: float = 0.2,
    eps_a: float = 0.05,
    eps_b: float = 0.05,
    I0: float = 0.0,
    A: float = 1.0,
    omega: float = 1.0,
) -> System:
    """Return a neuron that adapts its gain and threshold so its firing rate approaches an exponential distribution.

    Polyhomeostatic optimization targets the whole distribution of the firing rate, the one of maximal entropy for a
    given mean `mu`, rather than its mean alone. The neuron is driven by a periodic input, and the time average of its
    rate approaches `mu` as gain and threshold adapt.

    Args:
        mu:
            The target mean firing rate.

        eps_a:
            The adaptation rate of the gain.

        eps_b:
            The adaptation rate of the threshold.

        I0:
            The mean input.

        A:
            The amplitude of the periodic input.

        omega:
            The angular frequency of the periodic input.

    Returns:
        system:
            The system, with the potential `neuron.x`, the gain `neuron.a`, and the threshold `neuron.b`.

    Examples:
        >>> from skcomplex.systems.biology.neuroscience import polyhomeostatic_neuron
        >>> polyhomeostatic_neuron(mu=0.1).elements[0].mu
        0.1
    """
    neuron = PolyhomeostaticNeuronElement('neuron', mu=mu, eps_a=eps_a, eps_b=eps_b, I0=I0, A=A, omega=omega)
    interactions = [
        RelaxationInteraction(
            'membrane',
            target='neuron',
            state='x',
            goal=lambda s: s.target.I0 + s.target.A * sin(s.target.omega * s.t),
            time=lambda s: 1 / s.target.Gamma,
        ),
        RateInteraction(
            'gain_plasticity',
            target='neuron',
            state='a',
            rate=lambda s: s.target.eps_a * (1 / s.target.a + (s.target.x - s.target.b) * _kl_gradient(s)),
        ),
        RateInteraction(
            'threshold_plasticity',
            target='neuron',
            state='b',
            rate=lambda s: -s.target.eps_b * s.target.a * _kl_gradient(s),
        ),
    ]
    return System(elements=[neuron], interactions=interactions)

stdp_network(n_neurons=20, *, p=0.3, eta=0.5, spread=0.2, epsilon=0.05, A_plus=1.0, A_minus=1.05, tau=2.0, kappa=0.5, seed=0)

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.

Parameters:

Name	Type	Description	Default
n_neurons	int	The number of neurons.	20
p	float	The probability of every directed synapse.	0.3
eta	float	The mean excitability.	0.5
spread	float	The half-width of the uniform spread of the excitabilities.	0.2
epsilon	float	The learning rate.	0.05
A_plus	float	The amplitude of potentiation.	1.0
A_minus	float	The amplitude of depression.	1.05
tau	float	The time of both windows of the plasticity.	2.0
kappa	float	The initial efficacy of every synapse.	0.5
seed	int | None	The seed of the network, the excitabilities, and the initial phases.	0

Returns:

Name	Type	Description
system	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

Source code in src/skcomplex/systems/biology/neuroscience/_plastic_networks.py

@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])

theta_neurons(n_neurons=50, *, eta=-0.2, delta=0.1, kappa=1.5, n=2, random_state=0)

Return a network of theta neurons, equivalently quadratic integrate-and-fire neurons, coupled all to all.

Neuron i follows theta_i' = (1 - cos theta_i) + (1 + cos theta_i)(eta_i + kappa I), where I is the mean pulse a_n (1 - cos theta_j)^n of all neurons at the current state. The excitabilities are evenly spaced quantiles of a Lorentzian distribution of center eta and half-width delta, the distribution for which Ott and Antonsen's ansatz reduces the network exactly to one complex equation. Strong excitatory coupling switches the network from a quiet state to collective firing.

Parameters:

Name	Type	Description	Default
n_neurons	int	The number of neurons.	50
eta	float	The center of the excitabilities.	-0.2
delta	float	The half-width of the excitabilities.	0.1
kappa	float	The coupling strength.	1.5
n	int	The sharpness of the pulses, a positive integer.	2
random_state	int | Generator | None	The seed or generator of the initial phases.	0

Returns:

Name	Type	Description
system	System	The system, with the phases neurons.theta as its only state.

Examples:

>>> from skcomplex.systems.biology.neuroscience import theta_neurons
>>> theta_neurons(20).elements[0].prepare().size_
20

Source code in src/skcomplex/systems/biology/neuroscience/_populations.py

@register_system(
    field='biology/neuroscience',
    tags=('neuron', 'phase oscillators', 'mean field', 'network'),
    references=('Ermentrout and Kopell 1986', 'Luke, Barreto, and So 2013', 'arXiv:1902.05307'),
    defaults={'until': 20.0, 'dt': 0.1},
)
def theta_neurons(
    n_neurons: int = 50,
    *,
    eta: float = -0.2,
    delta: float = 0.1,
    kappa: float = 1.5,
    n: int = 2,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return a network of theta neurons, equivalently quadratic integrate-and-fire neurons, coupled all to all.

    Neuron `i` follows `theta_i' = (1 - cos theta_i) + (1 + cos theta_i)(eta_i + kappa I)`, where `I` is the mean pulse
    `a_n (1 - cos theta_j)^n` of all neurons at the current state. The excitabilities are evenly spaced quantiles of a
    Lorentzian distribution of center `eta` and half-width `delta`, the distribution for which Ott and Antonsen's
    ansatz reduces the network exactly to one complex equation. Strong excitatory coupling switches the network from a
    quiet state to collective firing.

    Args:
        n_neurons:
            The number of neurons.

        eta:
            The center of the excitabilities.

        delta:
            The half-width of the excitabilities.

        kappa:
            The coupling strength.

        n:
            The sharpness of the pulses, a positive integer.

        random_state:
            The seed or generator of the initial phases.

    Returns:
        system:
            The system, with the phases `neurons.theta` as its only state.

    Examples:
        >>> from skcomplex.systems.biology.neuroscience import theta_neurons
        >>> theta_neurons(20).elements[0].prepare().size_
        20
    """
    phases = np.random.default_rng(random_state).uniform(0.0, 2 * math.pi, n_neurons)
    neurons = GroupElement(
        'neurons',
        element=ThetaNeuronElement,
        size=n_neurons,
        values={'theta': phases, 'eta': np.asarray(_lorentzian_quantiles(eta, delta, n_neurons))},
    )
    return System(
        elements=[neurons],
        interactions=[PulseCouplingInteraction('coupling', neurons='neurons', kappa=kappa, n=n)],
    )

wilson_cowan(E=0.1, I=0.05, *, P=1.25, Q=0.0, c1=16.0, c2=12.0, c3=15.0, c4=3.0)

Return the Wilson-Cowan system of interacting excitatory and inhibitory populations of neurons.

With the published couplings, a moderate excitatory input P = 1.25 gives a limit cycle of population activity, and no input lets the activity decay to rest.

Parameters:

Name	Type	Description	Default
E	float	The initial excitatory activity.	0.1
I	float	The initial inhibitory activity.	0.05
P	float	The external input to the excitatory population.	1.25
Q	float	The external input to the inhibitory population.	0.0
c1	float	The excitatory to excitatory coupling.	16.0
c2	float	The inhibitory to excitatory coupling.	12.0
c3	float	The excitatory to inhibitory coupling.	15.0
c4	float	The inhibitory to inhibitory coupling.	3.0

Returns:

Name	Type	Description
system	System	The system, with the activities cortex.E and cortex.I.

Examples:

>>> from skcomplex.systems.biology.neuroscience import wilson_cowan
>>> wilson_cowan(P=0.0).elements[0].P
0.0

Source code in src/skcomplex/systems/biology/neuroscience/_populations.py

@register_system(
    field='biology/neuroscience',
    tags=('neural mass', 'oscillations', 'excitation-inhibition'),
    references=('Wilson and Cowan 1972', 'arXiv:2112.12147'),
    defaults={'until': 100.0, 'dt': 0.2},
)
def wilson_cowan(
    E: float = 0.1,
    I: float = 0.05,
    *,
    P: float = 1.25,
    Q: float = 0.0,
    c1: float = 16.0,
    c2: float = 12.0,
    c3: float = 15.0,
    c4: float = 3.0,
) -> System:
    """Return the Wilson-Cowan system of interacting excitatory and inhibitory populations of neurons.

    With the published couplings, a moderate excitatory input `P = 1.25` gives a limit cycle of population activity,
    and no input lets the activity decay to rest.

    Args:
        E:
            The initial excitatory activity.

        I:
            The initial inhibitory activity.

        P:
            The external input to the excitatory population.

        Q:
            The external input to the inhibitory population.

        c1:
            The excitatory to excitatory coupling.

        c2:
            The inhibitory to excitatory coupling.

        c3:
            The excitatory to inhibitory coupling.

        c4:
            The inhibitory to inhibitory coupling.

    Returns:
        system:
            The system, with the activities `cortex.E` and `cortex.I`.

    Examples:
        >>> from skcomplex.systems.biology.neuroscience import wilson_cowan
        >>> wilson_cowan(P=0.0).elements[0].P
        0.0
    """
    cortex = WilsonCowanElement('cortex', E=E, I=I, P=P, Q=Q, c1=c1, c2=c2, c3=c3, c4=c4)
    relaxations = [
        RelaxationInteraction(
            'excitatory_relaxation',
            target='cortex',
            state='E',
            goal=_excitation,
            time=lambda s: s.target.tau_e,
        ),
        RelaxationInteraction(
            'inhibitory_relaxation',
            target='cortex',
            state='I',
            goal=_inhibition,
            time=lambda s: s.target.tau_i,
        ),
    ]
    return System(elements=[cortex], interactions=relaxations)

winners_take_all(inputs=(1.0, 1.1, 0.9, 1.05, 0.95), *, beta=2.0, tau=1.0, x=0.1)

Return a winner-take-all network of units with self-limited growth and mutual inhibition.

Unit i follows tau x_i' = x_i (I_i - x_i - beta sum_{j != i} x_j), the competitive network of Fukai and Tanaka. One inhibition interaction binds every unit as an element set and reads the total activity of the set, so the network holds one interaction whatever the number of units. When the inhibition beta exceeds 1, only the unit with the largest input stays active, at activity I_max. When it is below 1, several winners share the activity. The winners-take-all network of Gros (2015, section 10.4), with reservoir variables and clique encoding, is a different model.

Parameters:

Name	Type	Description	Default
inputs	Sequence[float]	The input of every unit.	(1.0, 1.1, 0.9, 1.05, 0.95)
beta	float	The strength of the mutual inhibition.	2.0
tau	float	The time constant of the units.	1.0
x	float	The initial activity of every unit.	0.1

Returns:

Name	Type	Description
system	System	The system, with the activities unit_<i>.x.

Examples:

>>> from skcomplex.systems.biology.neuroscience import winners_take_all
>>> len(winners_take_all((1.0, 2.0, 3.0)).interactions)
1

Source code in src/skcomplex/systems/biology/neuroscience/_populations.py

@register_system(
    field='biology/neuroscience',
    tags=('neural network', 'competition', 'winner-take-all'),
    references=('Fukai and Tanaka 1997',),
    defaults={'until': 40.0, 'dt': 0.2},
)
def winners_take_all(
    inputs: Sequence[float] = (1.0, 1.1, 0.9, 1.05, 0.95),
    *,
    beta: float = 2.0,
    tau: float = 1.0,
    x: float = 0.1,
) -> System:
    """Return a winner-take-all network of units with self-limited growth and mutual inhibition.

    Unit `i` follows `tau x_i' = x_i (I_i - x_i - beta sum_{j != i} x_j)`, the competitive network of Fukai and
    Tanaka. One inhibition interaction binds every unit as an element set and reads the total activity of the set, so
    the network holds one interaction whatever the number of units. When the inhibition `beta` exceeds 1, only the unit
    with the largest input stays active, at activity `I_max`. When it is below 1, several winners share the activity.
    The winners-take-all network of Gros (2015, section 10.4), with reservoir variables and clique encoding, is a
    different model.

    Args:
        inputs:
            The input of every unit.

        beta:
            The strength of the mutual inhibition.

        tau:
            The time constant of the units.

        x:
            The initial activity of every unit.

    Returns:
        system:
            The system, with the activities `unit_<i>.x`.

    Examples:
        >>> from skcomplex.systems.biology.neuroscience import winners_take_all
        >>> len(winners_take_all((1.0, 2.0, 3.0)).interactions)
        1
    """
    units = [RateUnitElement(f'unit_{index}', x=x, I=value, tau=tau) for index, value in enumerate(inputs)]
    inhibition = MutualInhibitionInteraction('inhibition', units=[unit.label for unit in units], beta=beta)
    return System(elements=units, interactions=[inhibition])