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.models.biology.neuroscience import AdaptrodeElement
>>> AdaptrodeElement('synapse').w_max
1.0

Source code in src/skcomplex/models/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

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.models.biology.neuroscience import DelayedSynapseInteraction
>>> DelayedSynapseInteraction('synapse', source='a', target='b', J=-1.0).D
1.5

Source code in src/skcomplex/models/biology/neuroscience/_delayed_synapse.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.models.biology.neuroscience import ExcitableElement
>>> ExcitableElement('neuron').p
0.1

Source code in src/skcomplex/models/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.models.biology.neuroscience import FitzHughNagumoNeuronElement
>>> FitzHughNagumoNeuronElement('neuron').w
-0.6

Source code in src/skcomplex/models/biology/neuroscience/_fitz_hugh_nagumo_neuron.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.models.biology.neuroscience import GlobalInhibitionInteraction
>>> GlobalInhibitionInteraction('inhibition', inhibitor='z', oscillators='grid').W_z
1.5

Source code in src/skcomplex/models/biology/neuroscience/_global_inhibition.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.models.biology.neuroscience import HodgkinHuxleyElement
>>> HodgkinHuxleyElement('axon', I=10.0).g_na
120.0

Source code in src/skcomplex/models/biology/neuroscience/_hodgkin_huxley.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.models.biology.neuroscience import InhibitorDriveInteraction
>>> InhibitorDriveInteraction('drive', inhibitor='z', oscillators='grid').reduce
'max'

Source code in src/skcomplex/models/biology/neuroscience/_inhibitor_drive.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.models.biology.neuroscience import InhibitorElement
>>> InhibitorElement('inhibitor').phi
3.0

Source code in src/skcomplex/models/biology/neuroscience/_inhibitor.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.models.biology.neuroscience import IzhikevichElement
>>> IzhikevichElement('neuron').b
0.2

Source code in src/skcomplex/models/biology/neuroscience/_izhikevich.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.models.biology.neuroscience import LeakyIntegrateAndFireElement
>>> LeakyIntegrateAndFireElement('neuron').tau
10.0

Source code in src/skcomplex/models/biology/neuroscience/_leaky_integrate_and_fire.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.models.biology.neuroscience import MutualInhibitionInteraction
>>> MutualInhibitionInteraction('inhibition', units=['a', 'b'], beta=2.0).beta
2.0

Source code in src/skcomplex/models/biology/neuroscience/_mutual_inhibition.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.models.biology.neuroscience import NeuralFieldSiteElement
>>> NeuralFieldSiteElement('site').h
-0.2

Source code in src/skcomplex/models/biology/neuroscience/_neural_field_site.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

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.models.biology.neuroscience import PolyhomeostaticNeuronElement
>>> PolyhomeostaticNeuronElement('neuron').mu
0.2

Source code in src/skcomplex/models/biology/neuroscience/_polyhomeostatic_neuron.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.models.biology.neuroscience import PulseCouplingInteraction
>>> PulseCouplingInteraction('coupling', neurons='neurons', kappa=2.0).kappa
2.0

Source code in src/skcomplex/models/biology/neuroscience/_pulse_coupling.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.models.biology.neuroscience import RateUnitElement
>>> RateUnitElement('unit', I=1.2).I
1.2

Source code in src/skcomplex/models/biology/neuroscience/_rate_unit.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.elements import CompartmentElement
>>> from skcomplex.core.interactions import RateInteraction
>>> from skcomplex.core.system import System
>>> from skcomplex.engine.runtime import Simulation
>>> from skcomplex.models.biology.neuroscience import SpikeResetInteraction
>>> drive = RateInteraction('drive', target='neuron', state='n', rate=1.0)
>>> spike = SpikeResetInteraction('spike', target='neuron', state='n', threshold=1.0, after=0.0)
>>> system = System(elements=[CompartmentElement('neuron')], interactions=[drive, spike])
>>> trajectory = Simulation(system, dt=0.25).check().prepare().run(until=2.1).trajectory_
>>> [round(event.time, 6) for event in trajectory.events]
[1.0, 2.0]

Source code in src/skcomplex/models/biology/neuroscience/_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

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.models.biology.neuroscience import SpikingNeuronElement
>>> SpikingNeuronElement('neuron').tau_s
2.0

Source code in src/skcomplex/models/biology/neuroscience/_spiking_neuron.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

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.models.biology.neuroscience import TermanWangElement
>>> TermanWangElement('oscillator').gamma
6.0

Source code in src/skcomplex/models/biology/neuroscience/_terman_wang.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

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.elements import GroupElement
>>> from skcomplex.models.biology.neuroscience import ThetaNeuronElement
>>> GroupElement('neurons', element=ThetaNeuronElement, size=10, values={'eta': 0.5}).size
10

Source code in src/skcomplex/models/biology/neuroscience/_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

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.models.biology.neuroscience import WilsonCowanElement
>>> WilsonCowanElement('cortex', P=1.25).c1
16.0

Source code in src/skcomplex/models/biology/neuroscience/_wilson_cowan.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

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 model, with the weights synapse.w0, synapse.w1, and synapse.w2.

Examples:

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

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

@register_model(
    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 model, with the weights `synapse.w0`, `synapse.w1`, and `synapse.w2`.

    Examples:
        >>> from skcomplex.models.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 model, with the potential field.u.

Examples:

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

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

@register_model(
    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 model, with the potential `field.u`.

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

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

Examples:

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

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

@register_model(
    field='biology/neuroscience',
    tags=('neuron', 'spiking', 'network', 'delay', 'events'),
    references=('Brunel and Hakim 1999', 'arXiv:2112.12147'),
    defaults={'until': 100.0, 'dt': 0.5},
)
def delayed_lif_network(
    n_neurons: int = 10,
    *,
    probability: float = 0.3,
    J: float = 0.5,
    D: float = 1.5,
    mu: float = 25.0,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return a random network of leaky integrate-and-fire neurons whose spikes reach their targets after a delay.

    Every neuron is driven above threshold and fires on its own with period `tau ln((mu - reset) / (mu -
    threshold))`. Each directed link exists with a probability and moves the target's potential by `J` per spike,
    after the delay `D`. Excitation speeds the firing up and inhibition slows it down.

    Args:
        n_neurons:
            The number of neurons.

        probability:
            The probability of each directed link.

        J:
            The efficacy of every link in mV.

        D:
            The transmission delay in ms.

        mu:
            The constant drive in mV.

        random_state:
            The seed or generator of the links and the initial potentials.

    Returns:
        system:
            The model, with the potentials `neuron_<i>.v`, the traces `neuron_<i>.s`, and the spikes of neuron `i` as
            `spike_<i>` events.

    Examples:
        >>> from skcomplex.models.biology.neuroscience import delayed_lif_network
        >>> len(delayed_lif_network(5).elements)
        5
    """
    rng = np.random.default_rng(random_state)
    potentials = rng.uniform(10.0, 20.0, n_neurons)
    neurons = [
        SpikingNeuronElement(f'neuron_{index}', v=float(potentials[index]), mu=mu, tau_s=TRACE_TIME)
        for index in range(n_neurons)
    ]
    leaks = [
        RelaxationInteraction(
            f'leak_{index}',
            target=neuron.label,
            state='v',
            goal=lambda s: s.target.mu,
            time=lambda s: s.target.tau,
        )
        for index, neuron in enumerate(neurons)
    ]
    traces = [
        RateInteraction(
            f'trace_decay_{index}',
            target=neuron.label,
            state='s',
            rate=lambda s: -s.target.s / s.target.tau_s,
        )
        for index, neuron in enumerate(neurons)
    ]
    spikes = [
        SpikeResetInteraction(
            f'spike_{index}',
            target=neuron.label,
            threshold=THRESHOLD,
            after=RESET,
            kicked='s',
            kick=1 / TRACE_TIME,
        )
        for index, neuron in enumerate(neurons)
    ]
    synapses = [
        DelayedSynapseInteraction(
            f'synapse_{source}_{target}',
            source=f'neuron_{source}',
            target=f'neuron_{target}',
            J=J,
            D=D,
        )
        for source, target in _random_connections(n_neurons, probability, rng)
    ]
    return System(elements=neurons, interactions=[*leaks, *traces, *spikes, *synapses])

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 model, with the potential neuron.v and the recovery neuron.w.

Examples:

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

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

@register_model(
    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 model, with the potential `neuron.v` and the recovery `neuron.w`.

    Examples:
        >>> from skcomplex.models.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 model, with the potential neuron.v, the gates neuron.n, neuron.m, and neuron.h, and its spikes as spike events.

Examples:

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

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

@register_model(
    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 model, with the potential `neuron.v`, the gates `neuron.n`, `neuron.m`, and `neuron.h`, and its spikes
            as `spike` events.

    Examples:
        >>> from skcomplex.models.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 model, with the membrane potential neuron.v and its spikes recorded as spike events.

Examples:

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

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

@register_model(
    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 model, with the membrane potential `neuron.v` and its spikes recorded as `spike` events.

    Examples:
        >>> from skcomplex.models.biology.neuroscience import integrate_and_fire
        >>> from skcomplex.engine.runtime import Simulation
        >>> trajectory = Simulation(integrate_and_fire(), dt=0.1).check().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 model, with the potential neuron.v, the recovery neuron.u, and its spikes as spike events.

Examples:

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

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

@register_model(
    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 model, with the potential `neuron.v`, the recovery `neuron.u`, and its spikes as `spike` events.

    Examples:
        >>> from skcomplex.models.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 model, with the states neurons.phase, one value per node. Its standard observable activity is the fraction of excited elements, the firing rate per step.

Examples:

>>> import networkx as nx
>>> from skcomplex.models.biology.neuroscience import kinouchi_copelli
>>> kinouchi_copelli(nx.random_regular_graph(4, 20, seed=0), sigma=1.0).elements[0].values['p']
0.25

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

@register_model(
    field='biology/neuroscience',
    tags=('excitable', 'network', 'criticality', 'dynamic range', 'stochastic'),
    references=('Kinouchi and Copelli 2006', 'Larremore et al. 2011'),
    defaults={'until': 1000, 'dt': 1},
    observables=lambda system: [
        ObservableSpec('activity', group='neurons', variable='phase', reduce=_excited_fraction),
    ],
)
def kinouchi_copelli(
    graph: nx.Graph | None = None,
    *,
    sigma: float = 1.0,
    stimulus: float = 0.0,
    n_states: int = 5,
    excited: float | None = None,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return the Kinouchi-Copelli network of excitable elements, whose dynamic range peaks at the critical point.

    Every element excites each quiescent neighbour with the transmission probability `p = sigma / K`, where `K` is
    the mean degree of the graph, so `sigma` is the branching ratio, the mean number of elements one excitation
    excites. Below `sigma = 1` activity dies out without a stimulus, and above it activity sustains itself. The
    response, the mean firing rate against the stimulus rate, spans the most decades of stimulus at `sigma = 1`.

    Args:
        graph:
            The NetworkX graph, or `None` for an Erdős-Rényi graph of 1000 nodes with mean degree 10.

        sigma:
            The branching ratio.

        stimulus:
            The rate `r` of the external stimulus per step, which excites a quiescent element with probability
            `1 - exp(-r)`.

        n_states:
            The number of states of the cycle, one quiescent, one excited, and `n_states - 2` refractory.

        excited:
            The initial fraction of excited elements, with the others quiescent, or `None` for states drawn uniformly
            over the cycle.

        random_state:
            The seed or generator of the default graph and of the initial states.

    Returns:
        system:
            The model, with the states `neurons.phase`, one value per node. Its standard observable `activity` is the
            fraction of excited elements, the firing rate per step.

    Examples:
        >>> import networkx as nx
        >>> from skcomplex.models.biology.neuroscience import kinouchi_copelli
        >>> kinouchi_copelli(nx.random_regular_graph(4, 20, seed=0), sigma=1.0).elements[0].values['p']
        0.25
    """
    rng = np.random.default_rng(random_state)
    graph = nx.gnp_random_graph(1000, 10 / 999, seed=int(rng.integers(2**31))) if graph is None else graph
    size = graph.number_of_nodes()
    degree = 2 * graph.number_of_edges() / size
    phase = (
        rng.integers(0, n_states, size=size).astype(float)
        if excited is None
        else (rng.random(size) < excited).astype(float)
    )
    neurons = LatticeElement(
        'neurons',
        element=ExcitableElement,
        space=NetworkSpace(graph),
        values={'phase': phase, 'p': sigma / degree, 'stimulus': stimulus},
    )
    excitation = TransitionInteraction(
        'excitation',
        target='neurons',
        state='phase',
        source=0,
        destination=EXCITED,
        probability=_get_excitation,
    )
    return System(elements=[neurons], interactions=[excitation, *_get_cycle(n_states)], time=DiscreteTime())

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 model, with the fast variables oscillators.x, the slow variables oscillators.y, and the inhibitor inhibitor.z.

Examples:

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

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

@register_model(
    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 model, with the fast variables `oscillators.x`, the slow variables `oscillators.y`, and the inhibitor
            `inhibitor.z`.

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

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 model, with the potential neuron.x, the gain neuron.a, and the threshold neuron.b.

Examples:

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

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

@register_model(
    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 model, with the potential `neuron.x`, the gain `neuron.a`, and the threshold `neuron.b`.

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

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 model, with the phases neurons.theta as its only state.

Examples:

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

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

@register_model(
    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 model, with the phases `neurons.theta` as its only state.

    Examples:
        >>> from skcomplex.models.biology.neuroscience import theta_neurons
        >>> theta_neurons(20).elements[0].check().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 model 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 model, with the activities cortex.E and cortex.I.

Examples:

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

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

@register_model(
    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 model 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 model, with the activities `cortex.E` and `cortex.I`.

    Examples:
        >>> from skcomplex.models.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.

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 model, with the activities unit_<i>.x.

Examples:

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

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

@register_model(
    field='biology/neuroscience',
    tags=('neural network', 'competition', 'winner-take-all'),
    references=('Fukai and Tanaka 1997', 'Gros 2015'),
    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.

    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 model, with the activities `unit_<i>.x`.

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