Skip to content

Spatial triadic percolation on a torus

Spatial triadic percolation, where nearby signed regulators make the giant component form moving patterns.

  • Name: spatial_triadic_percolation
  • Field: mathematics/networks
  • Formalism: discrete map
  • References: Millán, Sun, Torres, and Bianconi 2023; Waxman 1988

Select this system

from skcomplex.systems import create_simulation, create_system, describe

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

Factory

@register_system(
    field='mathematics/networks',
    tags=('network', 'hypergraph', 'percolation', 'triadic interactions', 'spatial network', 'pattern formation'),
    references=('Millán, Sun, Torres, and Bianconi 2023', 'Waxman 1988'),
    defaults={'until': 100, 'dt': 1},
    observables=_giant_observables,
)
def spatial_triadic_percolation(
    side: float = 4.0,
    *,
    density: float = 100.0,
    c: float = 0.4,
    positive: float = 0.2,
    negative: float = 0.2,
    structural_range: float = 0.25,
    regulatory_range: float = 0.25,
    p: float = 0.8,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return spatial triadic percolation, where nearby signed regulators make the giant component form moving patterns.

    The nodes lie uniformly on a square torus of side `side` with `density` nodes per unit area, and two nodes at the
    distance `d` are linked with probability `c exp(-d / structural_range)`, a Waxman network, Eq. 1 of Millán et al.
    A node regulates a link positively with probability `positive exp(-d / regulatory_range)` and negatively with
    probability `negative exp(-d / regulatory_range)`, never both, from its distance `d` to the midpoint of the link,
    Eq. 2. Every step of the paper runs in two phases. In the phase `regulation` a link is active with probability `p`
    when none of its negative regulators is active and one of its positive regulators is, or, for the rare link
    without positive regulators, when none of its negative regulators is active. In the phase `percolation` a node is
    active when it belongs to the giant component of the active links. All nodes and links start active.

    The regulation is local, so the giant component forms clusters, stripes, and octopus-like patterns that move,
    and its size keeps changing instead of settling. At `p = 0.8` with `positive = negative` the patterns blink: a
    pattern deactivates the links it regulates negatively and the next step activates the region around it, so the
    active nodes of consecutive steps are anti-correlated. For long ranges the model becomes the triadic percolation of
    Erdős-Rényi networks with many regulators per link.

    Args:
        side:
            The side of the square torus.

        density:
            The number of nodes per unit area.

        c:
            The prefactor of the probability of a structural link.

        positive:
            The prefactor of the probability of a positive regulation.

        negative:
            The prefactor of the probability of a negative regulation.

        structural_range:
            The typical length of the structural links, `d0`.

        regulatory_range:
            The typical length of the regulations, `dr`.

        p:
            The probability that a link the regulation allows is active.

        random_state:
            The seed or generator of the positions, the links, and the regulators.

    Returns:
        system:
            The system, with the activities `nodes.active` and `links.active`, and `links.negatives` and
            `links.positives`, the numbers of negative and positive regulators of every link, on a network whose nodes
            hold their positions.

    Examples:
        >>> from skcomplex.core.engine.runtime import Simulation
        >>> from skcomplex.systems.mathematics.networks import spatial_triadic_percolation
        >>> system = spatial_triadic_percolation(1.0)
        >>> trajectory = Simulation(system, random_state=0).prepare().run(until=2).trajectory_
        >>> giant = trajectory.values('nodes.active').mean(axis=1)
        >>> bool(giant[0] == 1.0 and giant[1] < 0.1)
        True
    """
    rng = np.random.default_rng(random_state)
    graph = _waxman_torus(round(density * side**2), side, c, structural_range, rng)
    hyperedges, negatives, positives = _spatial_regulators(graph, side, (positive, negative), regulatory_range, rng)
    space = NetworkSpace(graph, hyperedges=hyperedges)
    nodes = LatticeElement('nodes', element=PercolatingNodeElement, space=space)
    values = {'negatives': negatives, 'positives': positives}
    links = GroupElement('links', element=SignedLinkElement, space=space, on='edges', values=values)
    interactions = [
        SignedRegulationInteraction('regulation', nodes='nodes', links='links', p=p, phase='regulation'),
        PercolationInteraction('percolation', nodes='nodes', links='links', phase='percolation'),
    ]
    time = DiscreteTime(schedule=PhasesScheduleSpec((PhaseSpec('regulation'), PhaseSpec('percolation'))))
    return System(elements=[nodes, links], interactions=interactions, time=time)

Components

Elements:

  • PercolatingNodeElement, from Element
  • SignedLinkElement, from Element

Interactions:

  • SignedRegulationInteraction, from HyperedgeInteraction
  • PercolationInteraction, from Interaction
class PercolatingNodeElement(Element):
    """Node of a network with spatial triadic regulation, active when it belongs to the giant component.

    Args:
        label:
            The label of the node group.

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

        active:
            1 when the node is in the giant component of the active links, and 0 otherwise.

    Examples:
        >>> from skcomplex.systems.mathematics.networks import PercolatingNodeElement
        >>> PercolatingNodeElement('nodes').active
        1.0
    """

    def __init__(
        self,
        label: str,
        visual: VisualSpec | None = ACTIVITY_VISUAL,
        *,
        active: Annotated[float, StateSpec(domain=FiniteEnumeratedSet([0, 1]))] = 1.0,
    ) -> None:
        super().__init__(label, visual=visual)
        self.active = active

class SignedLinkElement(Element):
    """Link regulated by signed nodes, with its activity and the numbers of its negative and positive regulators.

    Args:
        label:
            The label of the link group.

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

        active:
            1 when the link is active, and 0 otherwise.

        negatives:
            The number of regulating nodes that let the link be active only while they are inactive, the first
            regulators of its hyperedge.

        positives:
            The number of regulating nodes that let the link be active while they are active, the regulators of its
            hyperedge after the negative ones.

    Examples:
        >>> from skcomplex.systems.mathematics.networks import SignedLinkElement
        >>> SignedLinkElement('links', negatives=2).negatives
        2
    """

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

class SignedRegulationInteraction(HyperedgeInteraction):
    """Regulation of every link by signed regulating nodes, which let it be active with probability `p` in its phase.

    Every hyperedge `(i, j, r_1, ..., r_m)` holds the two ends of the link `(i, j)` and its regulating nodes, first its
    `negatives` negative regulators and then its `positives` positive ones, which the link reads. Further members pad
    the hyperedges to one size and regulate nothing, so links with several numbers of regulators share one size and one
    law writes every link once. In every step the link becomes active with probability `p` when none of its negative
    regulators is active and at least one of its positive regulators is active, or, for a link without positive
    regulators, when none of its negative regulators is active, and inactive otherwise.

    Args:
        label:
            The label of the interaction, unique among the parts of its system.

        nodes:
            The path of the group on the nodes, whose activity the regulators read.

        links:
            The path of the group on the edges, whose activity the regulation sets.

        p:
            The probability that a link the regulation allows is active.

        phase:
            The phase of the step the regulation runs in, or `None` for the default phase.

    Examples:
        >>> from skcomplex.systems.mathematics.networks import SignedRegulationInteraction
        >>> SignedRegulationInteraction('regulation', nodes='nodes', links='links').p
        0.5
    """

    def __init__(
        self,
        label: str,
        *,
        nodes: Annotated[str, RoleSpec()],
        links: Annotated[str, RoleSpec()],
        p: Annotated[float, ParameterSpec(domain=RealSet.closed(0, 1))] = 0.5,
        phase: str | None = None,
    ) -> None:
        super().__init__(label)
        self.nodes = nodes
        self.links = links
        self.p = p
        self.phase = phase

    def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
        """Return the activity of the link from its regulating nodes, the members of the hyperedge after its ends."""
        negatives, positives = s.edge.negatives, s.edge.positives
        blocked, silent = SR(1), SR(1)
        for position, regulator in enumerate(s.others[1:]):
            blocked *= 1 - where(negatives > position, regulator.active, 0)
            positive = where(negatives > position, 0, where(negatives + positives > position, 1, 0))
            silent *= 1 - positive * regulator.active
        allowed = blocked * where(positives > 0, 1 - silent, 1)
        return [UpdateLaw(s.edge.active, bernoulli('retained', s.p) * allowed, phase=self.phase)]

class PercolationInteraction(Interaction):
    """Activity of every node in its phase, 1 in the giant component of the active links and 0 outside it.

    Args:
        label:
            The label of the interaction, unique among the parts of its system.

        nodes:
            The path of the group on the nodes, whose activity the interaction sets.

        links:
            The path of the group on the edges, whose activity selects the links of the components.

        phase:
            The phase of the step the interaction runs in, or `None` for the default phase.

    Examples:
        >>> from skcomplex.systems.mathematics.networks import PercolationInteraction
        >>> PercolationInteraction('percolation', nodes='nodes', links='links').links
        'links'
    """

    def __init__(
        self,
        label: str,
        *,
        nodes: Annotated[str, RoleSpec()],
        links: Annotated[str, RoleSpec()],
        phase: str | None = None,
    ) -> None:
        super().__init__(label)
        self.nodes = nodes
        self.links = links
        self.phase = phase

    def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
        """Return the membership of the giant component of the active links."""
        return [UpdateLaw(s.nodes.active, components(s.links.active, reduce='giant'), phase=self.phase)]

Equations

nodes.active[i]: active_nodes(t + 1) == _neighbours_0
links.active[i]: active_links(t + 1) == (where(negatives_links__edge > 11, active_nodes__other12, 0) - 1)*(where(negatives_links__edge > 27, active_nodes__other28, 0) - 1)*(where(negatives_links__edge > 13, active_nodes__other14, 0) - 1)*(where(negatives_links__edge > 0, active_nodes__other1, 0) - 1)*(where(negatives_links__edge > 24, active_nodes__other25, 0) - 1)*(where(negatives_links__edge > 4, active_nodes__other5, 0) - 1)*(where(negatives_links__edge > 3, active_nodes__other4, 0) - 1)*(where(negatives_links__edge > 31, active_nodes__other32, 0) - 1)*(where(negatives_links__edge > 14, active_nodes__other15, 0) - 1)*(where(negatives_links__edge > 21, active_nodes__other22, 0) - 1)*(where(negatives_links__edge > 1, active_nodes__other2, 0) - 1)*(where(negatives_links__edge > 25, active_nodes__other26, 0) - 1)*(where(negatives_links__edge > 28, active_nodes__other29, 0) - 1)*(where(negatives_links__edge > 18, active_nodes__other19, 0) - 1)*(where(negatives_links__edge > 7, active_nodes__other8, 0) - 1)*(where(negatives_links__edge > 10, active_nodes__other11, 0) - 1)*(where(negatives_links__edge > 15, active_nodes__other16, 0) - 1)*(where(negatives_links__edge > 8, active_nodes__other9, 0) - 1)*(where(negatives_links__edge > 12, active_nodes__other13, 0) - 1)*(where(negatives_links__edge > 19, active_nodes__other20, 0) - 1)*(where(negatives_links__edge > 2, active_nodes__other3, 0) - 1)*(where(negatives_links__edge > 5, active_nodes__other6, 0) - 1)*(where(negatives_links__edge > 29, active_nodes__other30, 0) - 1)*(where(negatives_links__edge > 23, active_nodes__other24, 0) - 1)*(where(negatives_links__edge > 9, active_nodes__other10, 0) - 1)*(where(negatives_links__edge > 16, active_nodes__other17, 0) - 1)*(where(negatives_links__edge > 26, active_nodes__other27, 0) - 1)*(where(negatives_links__edge > 20, active_nodes__other21, 0) - 1)*(where(negatives_links__edge > 6, active_nodes__other7, 0) - 1)*(where(negatives_links__edge > 30, active_nodes__other31, 0) - 1)*(where(negatives_links__edge > 17, active_nodes__other18, 0) - 1)*(where(negatives_links__edge > 22, active_nodes__other23, 0) - 1)*where(positives_links__edge > 0, -(active_nodes__other20*where(negatives_links__edge > 19, 0, where(negatives_links__edge + positives_links__edge > 19, 1, 0)) - 1)*(active_nodes__other32*where(negatives_links__edge > 31, 0, where(negatives_links__edge + positives_links__edge > 31, 1, 0)) - 1)*(active_nodes__other24*where(negatives_links__edge > 23, 0, where(negatives_links__edge + positives_links__edge > 23, 1, 0)) - 1)*(active_nodes__other3*where(negatives_links__edge > 2, 0, where(negatives_links__edge + positives_links__edge > 2, 1, 0)) - 1)*(active_nodes__other6*where(negatives_links__edge > 5, 0, where(negatives_links__edge + positives_links__edge > 5, 1, 0)) - 1)*(active_nodes__other30*where(negatives_links__edge > 29, 0, where(negatives_links__edge + positives_links__edge > 29, 1, 0)) - 1)*(active_nodes__other13*where(negatives_links__edge > 12, 0, where(negatives_links__edge + positives_links__edge > 12, 1, 0)) - 1)*(active_nodes__other21*where(negatives_links__edge > 20, 0, where(negatives_links__edge + positives_links__edge > 20, 1, 0)) - 1)*(active_nodes__other16*where(negatives_links__edge > 15, 0, where(negatives_links__edge + positives_links__edge > 15, 1, 0)) - 1)*(active_nodes__other27*where(negatives_links__edge > 26, 0, where(negatives_links__edge + positives_links__edge > 26, 1, 0)) - 1)*(active_nodes__other17*where(negatives_links__edge > 16, 0, where(negatives_links__edge + positives_links__edge > 16, 1, 0)) - 1)*(active_nodes__other10*where(negatives_links__edge > 9, 0, where(negatives_links__edge + positives_links__edge > 9, 1, 0)) - 1)*(active_nodes__other18*where(negatives_links__edge > 17, 0, where(negatives_links__edge + positives_links__edge > 17, 1, 0)) - 1)*(active_nodes__other9*where(negatives_links__edge > 8, 0, where(negatives_links__edge + positives_links__edge > 8, 1, 0)) - 1)*(active_nodes__other7*where(negatives_links__edge > 6, 0, where(negatives_links__edge + positives_links__edge > 6, 1, 0)) - 1)*(active_nodes__other28*where(negatives_links__edge > 27, 0, where(negatives_links__edge + positives_links__edge > 27, 1, 0)) - 1)*(active_nodes__other19*where(negatives_links__edge > 18, 0, where(negatives_links__edge + positives_links__edge > 18, 1, 0)) - 1)*(active_nodes__other31*where(negatives_links__edge > 30, 0, where(negatives_links__edge + positives_links__edge > 30, 1, 0)) - 1)*(active_nodes__other14*where(negatives_links__edge > 13, 0, where(negatives_links__edge + positives_links__edge > 13, 1, 0)) - 1)*(active_nodes__other4*where(negatives_links__edge > 3, 0, where(negatives_links__edge + positives_links__edge > 3, 1, 0)) - 1)*(active_nodes__other15*where(negatives_links__edge > 14, 0, where(negatives_links__edge + positives_links__edge > 14, 1, 0)) - 1)*(active_nodes__other12*where(negatives_links__edge > 11, 0, where(negatives_links__edge + positives_links__edge > 11, 1, 0)) - 1)*(active_nodes__other5*where(negatives_links__edge > 4, 0, where(negatives_links__edge + positives_links__edge > 4, 1, 0)) - 1)*(active_nodes__other8*where(negatives_links__edge > 7, 0, where(negatives_links__edge + positives_links__edge > 7, 1, 0)) - 1)*(active_nodes__other25*where(negatives_links__edge > 24, 0, where(negatives_links__edge + positives_links__edge > 24, 1, 0)) - 1)*(active_nodes__other1*where(negatives_links__edge > 0, 0, where(negatives_links__edge + positives_links__edge > 0, 1, 0)) - 1)*(active_nodes__other22*where(negatives_links__edge > 21, 0, where(negatives_links__edge + positives_links__edge > 21, 1, 0)) - 1)*(active_nodes__other26*where(negatives_links__edge > 25, 0, where(negatives_links__edge + positives_links__edge > 25, 1, 0)) - 1)*(active_nodes__other2*where(negatives_links__edge > 1, 0, where(negatives_links__edge + positives_links__edge > 1, 1, 0)) - 1)*(active_nodes__other29*where(negatives_links__edge > 28, 0, where(negatives_links__edge + positives_links__edge > 28, 1, 0)) - 1)*(active_nodes__other23*where(negatives_links__edge > 22, 0, where(negatives_links__edge + positives_links__edge > 22, 1, 0)) - 1)*(active_nodes__other11*where(negatives_links__edge > 10, 0, where(negatives_links__edge + positives_links__edge > 10, 1, 0)) - 1) + 1, 1)*where(bernoulli_retained_regulation < p_regulation, 1, 0)