Theta neurons and QIF networks

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

• Name: theta_neurons
• Field: biology/neuroscience
• Formalism: ordinary differential equations
• References: Ermentrout and Kopell 1986; Luke, Barreto, and So 2013; arXiv:1902.05307

Factory

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

Components

class ThetaNeuronElement(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`.

    Args:
        label:
            The label of the neuron or the group.

        visual:
            How the element is drawn, or `None` for `default_visual_`.

        theta:
            The initial phase.

        eta:
            The excitability.

    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
    """

    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

    def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
        """Return the phase dynamics.

        Args:
            s:
                The namespace of the symbols the laws can use.

        Returns:
            laws:
                The laws of the component.
        """
        return [RateLaw(s.theta, 1 - cos(s.theta) + (1 + cos(s.theta)) * s.eta)]

class PulseCouplingInteraction(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.

    Args:
        label:
            The label of the interaction.

        neurons:
            The path of the group of theta neurons.

        kappa:
            The coupling strength.

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

    Examples:
        >>> from skcomplex.models.biology.neuroscience import PulseCouplingInteraction
        >>> PulseCouplingInteraction('coupling', neurons='neurons', kappa=2.0).kappa
        2.0
    """

    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

    def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
        """Return the current of the mean pulse into every neuron.

        Args:
            s:
                The namespace of the symbols the laws can use.

        Returns:
            laws:
                The laws of the component.
        """
        theta = s.neurons.theta
        normalization = 2**self.n * math.factorial(self.n) ** 2 / math.factorial(2 * self.n)
        pulses = mean(normalization * (1 - cos(theta)) ** self.n)
        return [RateLaw(theta, (1 + cos(theta)) * s.kappa * pulses)]

Equations

neurons.theta[i]: diff(theta_neurons(t), t) == _reduced_0*kappa_coupling*(cos(theta_neurons(t)) + 1) + eta_neurons*(cos(theta_neurons(t)) + 1) - cos(theta_neurons(t)) + 1