Skip to content

Phase units with activity-dependent conduction delays

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

  • Name: adaptive_delay_network
  • Field: biology/neuroscience
  • Formalism: ordinary differential equations
  • References: Park and Lefebvre 2020; Berner et al. 2023

Select this system

from skcomplex.systems import create_simulation, create_system, describe

system = create_system('adaptive_delay_network')
entry = describe('adaptive_delay_network')
trajectory = create_simulation('adaptive_delay_network').prepare().run().trajectory_

Factory

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

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

    Args:
        n_units:
            The number of units.

        omega:
            The mean natural frequency.

        spread:
            The standard deviation of the natural frequencies.

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

        epsilon:
            The adaptation rate of the delays.

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

        tau0:
            The baseline delay of every edge.

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

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

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

Components

Elements:

  • PhaseUnitElement, from Element
  • AxonElement, from Element

Interactions:

  • DelayAdaptationInteraction, from Interaction
class PhaseUnitElement(Element):
    """Neural unit reduced to a phase that turns at its natural frequency, `theta' = omega`.

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

    Args:
        label:
            The label of the unit or of the group of units.

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

        theta:
            The initial phase.

        omega:
            The natural frequency.

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

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

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

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

        Returns:
            laws:
                The laws of the component.
        """
        return [RateLaw(s.theta, s.omega)]

class AxonElement(Element):
    """Axon on an edge of a network, whose conduction delay myelination changes.

    Args:
        label:
            The label of the axon or of the group on the edges.

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

        tau:
            The initial conduction delay.

        tau0:
            The baseline delay that the delay relaxes to.

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

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

class DelayAdaptationInteraction(Interaction):
    """Phase coupling through conduction delays that adapt to the phase differences of the units they join.

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

    Args:
        label:
            The label of the interaction.

        source:
            The path of the unit that sends the signal.

        target:
            The path of the unit that receives the signal.

        axon:
            The path of the axon that conducts it.

        kappa:
            The coupling strength of an edge.

        epsilon:
            The adaptation rate of the delays.

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

        width:
            The width of the smooth step.

        shortest:
            The shortest delay the coupling reads.

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

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

    def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
        """Return the delayed coupling of the target and the adaptation of the delay.

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

        Returns:
            laws:
                The laws of the component.
        """
        source, target, tau = s.source.theta, s.target.theta, s.axon.tau
        step = (1 + tanh(tau / s.width)) / 2
        return [
            RateLaw(target, s.kappa * sin(delayed(source, clip(tau, low=self.shortest), self.shortest) - target)),
            RateLaw(tau, s.epsilon * step * (-(tau - s.axon.tau0) + s.rho * sin(source - target))),
        ]

Equations

unit_0.theta: diff(theta_unit_0(t), t) == kappa_myelination_1_0*sin(-theta_unit_0(t) + theta_unit_1(t - where(tau_axon_1_0(t) < 0.1, 0.1, tau_axon_1_0(t)))) + kappa_myelination_2_0*sin(-theta_unit_0(t) + theta_unit_2(t - where(tau_axon_2_0(t) < 0.1, 0.1, tau_axon_2_0(t)))) + kappa_myelination_3_0*sin(-theta_unit_0(t) + theta_unit_3(t - where(tau_axon_3_0(t) < 0.1, 0.1, tau_axon_3_0(t)))) + kappa_myelination_4_0*sin(-theta_unit_0(t) + theta_unit_4(t - where(tau_axon_4_0(t) < 0.1, 0.1, tau_axon_4_0(t)))) + omega_unit_0
unit_1.theta: diff(theta_unit_1(t), t) == kappa_myelination_0_1*sin(theta_unit_0(t - where(tau_axon_0_1(t) < 0.1, 0.1, tau_axon_0_1(t))) - theta_unit_1(t)) + kappa_myelination_2_1*sin(-theta_unit_1(t) + theta_unit_2(t - where(tau_axon_2_1(t) < 0.1, 0.1, tau_axon_2_1(t)))) + kappa_myelination_3_1*sin(-theta_unit_1(t) + theta_unit_3(t - where(tau_axon_3_1(t) < 0.1, 0.1, tau_axon_3_1(t)))) + kappa_myelination_4_1*sin(-theta_unit_1(t) + theta_unit_4(t - where(tau_axon_4_1(t) < 0.1, 0.1, tau_axon_4_1(t)))) + omega_unit_1
unit_2.theta: diff(theta_unit_2(t), t) == kappa_myelination_0_2*sin(theta_unit_0(t - where(tau_axon_0_2(t) < 0.1, 0.1, tau_axon_0_2(t))) - theta_unit_2(t)) + kappa_myelination_1_2*sin(theta_unit_1(t - where(tau_axon_1_2(t) < 0.1, 0.1, tau_axon_1_2(t))) - theta_unit_2(t)) + kappa_myelination_3_2*sin(-theta_unit_2(t) + theta_unit_3(t - where(tau_axon_3_2(t) < 0.1, 0.1, tau_axon_3_2(t)))) + kappa_myelination_4_2*sin(-theta_unit_2(t) + theta_unit_4(t - where(tau_axon_4_2(t) < 0.1, 0.1, tau_axon_4_2(t)))) + omega_unit_2
unit_3.theta: diff(theta_unit_3(t), t) == kappa_myelination_0_3*sin(theta_unit_0(t - where(tau_axon_0_3(t) < 0.1, 0.1, tau_axon_0_3(t))) - theta_unit_3(t)) + kappa_myelination_1_3*sin(theta_unit_1(t - where(tau_axon_1_3(t) < 0.1, 0.1, tau_axon_1_3(t))) - theta_unit_3(t)) + kappa_myelination_2_3*sin(theta_unit_2(t - where(tau_axon_2_3(t) < 0.1, 0.1, tau_axon_2_3(t))) - theta_unit_3(t)) + kappa_myelination_4_3*sin(-theta_unit_3(t) + theta_unit_4(t - where(tau_axon_4_3(t) < 0.1, 0.1, tau_axon_4_3(t)))) + omega_unit_3
unit_4.theta: diff(theta_unit_4(t), t) == kappa_myelination_0_4*sin(theta_unit_0(t - where(tau_axon_0_4(t) < 0.1, 0.1, tau_axon_0_4(t))) - theta_unit_4(t)) + kappa_myelination_1_4*sin(theta_unit_1(t - where(tau_axon_1_4(t) < 0.1, 0.1, tau_axon_1_4(t))) - theta_unit_4(t)) + kappa_myelination_2_4*sin(theta_unit_2(t - where(tau_axon_2_4(t) < 0.1, 0.1, tau_axon_2_4(t))) - theta_unit_4(t)) + kappa_myelination_3_4*sin(theta_unit_3(t - where(tau_axon_3_4(t) < 0.1, 0.1, tau_axon_3_4(t))) - theta_unit_4(t)) + omega_unit_4
axon_0_1.tau: diff(tau_axon_0_1(t), t) == 1/2*(rho_myelination_0_1*sin(theta_unit_0(t) - theta_unit_1(t)) + tau0_axon_0_1 - tau_axon_0_1(t))*epsilon_myelination_0_1*(tanh(tau_axon_0_1(t)/width_myelination_0_1) + 1)
axon_0_2.tau: diff(tau_axon_0_2(t), t) == 1/2*(rho_myelination_0_2*sin(theta_unit_0(t) - theta_unit_2(t)) + tau0_axon_0_2 - tau_axon_0_2(t))*epsilon_myelination_0_2*(tanh(tau_axon_0_2(t)/width_myelination_0_2) + 1)
axon_0_3.tau: diff(tau_axon_0_3(t), t) == 1/2*(rho_myelination_0_3*sin(theta_unit_0(t) - theta_unit_3(t)) + tau0_axon_0_3 - tau_axon_0_3(t))*epsilon_myelination_0_3*(tanh(tau_axon_0_3(t)/width_myelination_0_3) + 1)
axon_0_4.tau: diff(tau_axon_0_4(t), t) == 1/2*(rho_myelination_0_4*sin(theta_unit_0(t) - theta_unit_4(t)) + tau0_axon_0_4 - tau_axon_0_4(t))*epsilon_myelination_0_4*(tanh(tau_axon_0_4(t)/width_myelination_0_4) + 1)
axon_1_0.tau: diff(tau_axon_1_0(t), t) == 1/2*(rho_myelination_1_0*sin(-theta_unit_0(t) + theta_unit_1(t)) + tau0_axon_1_0 - tau_axon_1_0(t))*epsilon_myelination_1_0*(tanh(tau_axon_1_0(t)/width_myelination_1_0) + 1)
axon_1_2.tau: diff(tau_axon_1_2(t), t) == 1/2*(rho_myelination_1_2*sin(theta_unit_1(t) - theta_unit_2(t)) + tau0_axon_1_2 - tau_axon_1_2(t))*epsilon_myelination_1_2*(tanh(tau_axon_1_2(t)/width_myelination_1_2) + 1)
axon_1_3.tau: diff(tau_axon_1_3(t), t) == 1/2*(rho_myelination_1_3*sin(theta_unit_1(t) - theta_unit_3(t)) + tau0_axon_1_3 - tau_axon_1_3(t))*epsilon_myelination_1_3*(tanh(tau_axon_1_3(t)/width_myelination_1_3) + 1)
axon_1_4.tau: diff(tau_axon_1_4(t), t) == 1/2*(rho_myelination_1_4*sin(theta_unit_1(t) - theta_unit_4(t)) + tau0_axon_1_4 - tau_axon_1_4(t))*epsilon_myelination_1_4*(tanh(tau_axon_1_4(t)/width_myelination_1_4) + 1)
axon_2_0.tau: diff(tau_axon_2_0(t), t) == 1/2*(rho_myelination_2_0*sin(-theta_unit_0(t) + theta_unit_2(t)) + tau0_axon_2_0 - tau_axon_2_0(t))*epsilon_myelination_2_0*(tanh(tau_axon_2_0(t)/width_myelination_2_0) + 1)
axon_2_1.tau: diff(tau_axon_2_1(t), t) == 1/2*(rho_myelination_2_1*sin(-theta_unit_1(t) + theta_unit_2(t)) + tau0_axon_2_1 - tau_axon_2_1(t))*epsilon_myelination_2_1*(tanh(tau_axon_2_1(t)/width_myelination_2_1) + 1)
axon_2_3.tau: diff(tau_axon_2_3(t), t) == 1/2*(rho_myelination_2_3*sin(theta_unit_2(t) - theta_unit_3(t)) + tau0_axon_2_3 - tau_axon_2_3(t))*epsilon_myelination_2_3*(tanh(tau_axon_2_3(t)/width_myelination_2_3) + 1)
axon_2_4.tau: diff(tau_axon_2_4(t), t) == 1/2*(rho_myelination_2_4*sin(theta_unit_2(t) - theta_unit_4(t)) + tau0_axon_2_4 - tau_axon_2_4(t))*epsilon_myelination_2_4*(tanh(tau_axon_2_4(t)/width_myelination_2_4) + 1)
axon_3_0.tau: diff(tau_axon_3_0(t), t) == 1/2*(rho_myelination_3_0*sin(-theta_unit_0(t) + theta_unit_3(t)) + tau0_axon_3_0 - tau_axon_3_0(t))*epsilon_myelination_3_0*(tanh(tau_axon_3_0(t)/width_myelination_3_0) + 1)
axon_3_1.tau: diff(tau_axon_3_1(t), t) == 1/2*(rho_myelination_3_1*sin(-theta_unit_1(t) + theta_unit_3(t)) + tau0_axon_3_1 - tau_axon_3_1(t))*epsilon_myelination_3_1*(tanh(tau_axon_3_1(t)/width_myelination_3_1) + 1)
axon_3_2.tau: diff(tau_axon_3_2(t), t) == 1/2*(rho_myelination_3_2*sin(-theta_unit_2(t) + theta_unit_3(t)) + tau0_axon_3_2 - tau_axon_3_2(t))*epsilon_myelination_3_2*(tanh(tau_axon_3_2(t)/width_myelination_3_2) + 1)
axon_3_4.tau: diff(tau_axon_3_4(t), t) == 1/2*(rho_myelination_3_4*sin(theta_unit_3(t) - theta_unit_4(t)) + tau0_axon_3_4 - tau_axon_3_4(t))*epsilon_myelination_3_4*(tanh(tau_axon_3_4(t)/width_myelination_3_4) + 1)
axon_4_0.tau: diff(tau_axon_4_0(t), t) == 1/2*(rho_myelination_4_0*sin(-theta_unit_0(t) + theta_unit_4(t)) + tau0_axon_4_0 - tau_axon_4_0(t))*epsilon_myelination_4_0*(tanh(tau_axon_4_0(t)/width_myelination_4_0) + 1)
axon_4_1.tau: diff(tau_axon_4_1(t), t) == 1/2*(rho_myelination_4_1*sin(-theta_unit_1(t) + theta_unit_4(t)) + tau0_axon_4_1 - tau_axon_4_1(t))*epsilon_myelination_4_1*(tanh(tau_axon_4_1(t)/width_myelination_4_1) + 1)
axon_4_2.tau: diff(tau_axon_4_2(t), t) == 1/2*(rho_myelination_4_2*sin(-theta_unit_2(t) + theta_unit_4(t)) + tau0_axon_4_2 - tau_axon_4_2(t))*epsilon_myelination_4_2*(tanh(tau_axon_4_2(t)/width_myelination_4_2) + 1)
axon_4_3.tau: diff(tau_axon_4_3(t), t) == 1/2*(rho_myelination_4_3*sin(-theta_unit_3(t) + theta_unit_4(t)) + tau0_axon_4_3 - tau_axon_4_3(t))*epsilon_myelination_4_3*(tanh(tau_axon_4_3(t)/width_myelination_4_3) + 1)