Skip to content

Kauffman random Boolean networks

A Kauffman random Boolean network, nodes updated together by random functions of k random inputs.

  • Name: random_boolean_network
  • Field: biology/systems_biology
  • Formalism: discrete map
  • References: Kauffman 1969; Gros 2015

Factory

@register_model(
    field='biology/systems_biology',
    tags=('boolean network', 'gene regulation', 'discrete map', 'criticality'),
    references=('Kauffman 1969', 'Gros 2015'),
    defaults={'until': 50, 'dt': 1},
)
def random_boolean_network(
    n_nodes: int = 20,
    k: int = 2,
    *,
    bias: float = 0.5,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return a Kauffman random Boolean network, nodes updated together by random functions of `k` random inputs.

    Each node reads `k` distinct nodes, possibly itself, through a truth table whose entries are 1 with probability
    `bias`. The finite state space makes every run end on a cycle. With `k = 2` and `bias = 0.5` the network is
    critical, between frozen and chaotic.

    Args:
        n_nodes:
            The number of nodes.

        k:
            The number of inputs per node, at most the number of nodes.

        bias:
            The probability that an entry of a truth table is 1.

        random_state:
            The seed or generator of the wiring, the tables, and the initial values.

    Returns:
        system:
            The model, with the values `node_<i>.x`.

    Raises:
        ValidationError:
            `k` is not between 1 and the number of nodes.

    Examples:
        >>> from skcomplex.models.biology.systems_biology import random_boolean_network
        >>> len(random_boolean_network(10, 3).interactions[0].table)
        8
    """
    if not 1 <= k <= n_nodes:
        error_msg = f'k = {k} inputs per node must be between 1 and the number of nodes, n_nodes = {n_nodes}.'
        raise ValidationError(error_msg)
    rng = np.random.default_rng(random_state)
    nodes = [BooleanNodeElement(f'node_{index}', x=int(rng.integers(2))) for index in range(n_nodes)]
    functions = [
        BooleanFunctionInteraction(
            f'function_{index}',
            target=f'node_{index}',
            inputs=[f'node_{source}' for source in rng.choice(n_nodes, size=k, replace=False)],
            table=tuple(int(entry) for entry in rng.random(2**k) < bias),
        )
        for index in range(n_nodes)
    ]
    return System(elements=nodes, interactions=functions)

Components

class BooleanNodeElement(Element):
    """Node of a Boolean network, a gene or protein that is on or off.

    Args:
        label:
            The label of the node.

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

        x:
            The initial value, 0 or 1.

    Examples:
        >>> from skcomplex.models.biology.systems_biology import BooleanNodeElement
        >>> BooleanNodeElement('Sic1', x=1).x
        1
    """

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

class BooleanFunctionInteraction(MapInteraction):
    """Synchronous update of a node by a truth table of any number of input nodes.

    A table of `2^k` entries reads the `k` inputs as the bits of a binary number, from the most significant, and gives
    the next value of the target. The update is the sum, over the patterns whose entry is 1, of the product that is 1
    exactly on that pattern.

    Args:
        label:
            The label of the interaction.

        target:
            The path of the updated node.

        inputs:
            The paths of the input nodes, in the order of the bits of the table.

        table:
            The next value of the target for each of the `2^k` input patterns.

    Raises:
        ValidationError:
            The table does not have one entry per pattern of the inputs.

    Examples:
        >>> from skcomplex.models.biology.systems_biology import BooleanFunctionInteraction
        >>> BooleanFunctionInteraction('and', target='c', inputs=['a', 'b'], table=(0, 0, 0, 1)).table
        (0, 0, 0, 1)
    """

    def __init__(
        self,
        label: str,
        *,
        target: Annotated[str, RoleSpec()],
        inputs: Annotated[list[str], RoleSpec(variadic=True)],
        table: tuple[int, ...] = (0,),
    ) -> None:
        super().__init__(label, target=target, state='x', function=self._get_next)
        self.inputs = inputs
        self.table = table

    def _get_pattern(self, inputs: Sequence[Expression], index: int) -> Expression:
        """Return the product that is 1 when the inputs, read as bits from the most significant, spell an index."""
        bits = [(index >> (len(inputs) - 1 - position)) & 1 for position in range(len(inputs))]
        return prod((value if bit else 1 - value for value, bit in zip(inputs, bits, strict=True)), start=SR(1))

    def _get_next(self, s: SimpleNamespace) -> Expression:
        """Return the next value of the target by the truth table."""
        inputs = [item.x for item in s.inputs]
        terms = [self._get_pattern(inputs, index) for index, entry in enumerate(self.table) if entry]
        return sum(terms, start=SR(0))

    def _check(self) -> None:
        """Raise an error for a table without one entry per pattern of the inputs.

        Raises:
            ValidationError:
                The table does not have one entry per pattern of the inputs.
        """
        super()._check()
        if len(self.table) != 2 ** len(self.inputs):
            error_msg = (
                f'BooleanFunction {self.label} reads {len(self.inputs)} inputs, which take {2 ** len(self.inputs)} '
                f'patterns, but its table has {len(self.table)} entries.'
            )
            raise ValidationError(error_msg)

Equations

node_0.x: x_node_0(t + 1) == -(x_node_5(t) - 1)*x_node_16(t) + x_node_16(t)*x_node_5(t)
node_1.x: x_node_1(t + 1) == -(x_node_1(t) - 1)*x_node_0(t) - (x_node_0(t) - 1)*x_node_1(t) + x_node_0(t)*x_node_1(t)
node_2.x: x_node_2(t + 1) == -(x_node_0(t) - 1)*x_node_2(t)
node_3.x: x_node_3(t + 1) == x_node_12(t)*x_node_16(t)
node_4.x: x_node_4(t + 1) == -(x_node_2(t) - 1)*x_node_16(t) - (x_node_16(t) - 1)*x_node_2(t)
node_5.x: x_node_5(t + 1) == (x_node_1(t) - 1)*(x_node_13(t) - 1) - (x_node_1(t) - 1)*x_node_13(t)
node_6.x: x_node_6(t + 1) == -(x_node_9(t) - 1)*x_node_6(t) + x_node_6(t)*x_node_9(t)
node_7.x: x_node_7(t + 1) == -(x_node_7(t) - 1)*x_node_19(t) + x_node_19(t)*x_node_7(t)
node_8.x: x_node_8(t + 1) == (x_node_12(t) - 1)*(x_node_6(t) - 1) - (x_node_6(t) - 1)*x_node_12(t) + x_node_12(t)*x_node_6(t)
node_9.x: x_node_9(t + 1) == (x_node_11(t) - 1)*(x_node_2(t) - 1) - (x_node_11(t) - 1)*x_node_2(t) + x_node_11(t)*x_node_2(t)
node_10.x: x_node_10(t + 1) == -(x_node_17(t) - 1)*x_node_13(t) - (x_node_13(t) - 1)*x_node_17(t)
node_11.x: x_node_11(t + 1) == (x_node_15(t) - 1)*(x_node_3(t) - 1) - (x_node_15(t) - 1)*x_node_3(t) + x_node_15(t)*x_node_3(t)
node_12.x: x_node_12(t + 1) == -(x_node_12(t) - 1)*x_node_17(t) + x_node_12(t)*x_node_17(t)
node_13.x: x_node_13(t + 1) == -(x_node_18(t) - 1)*x_node_8(t)
node_14.x: x_node_14(t + 1) == x_node_17(t)*x_node_18(t)
node_15.x: x_node_15(t + 1) == -(x_node_19(t) - 1)*x_node_2(t)
node_16.x: x_node_16(t + 1) == (x_node_17(t) - 1)*(x_node_6(t) - 1) - (x_node_17(t) - 1)*x_node_6(t)
node_17.x: x_node_17(t + 1) == -(x_node_19(t) - 1)*x_node_15(t) + x_node_15(t)*x_node_19(t)
node_18.x: x_node_18(t + 1) == (x_node_14(t) - 1)*(x_node_15(t) - 1) - (x_node_15(t) - 1)*x_node_14(t)
node_19.x: x_node_19(t + 1) == -(x_node_11(t) - 1)*x_node_2(t) + x_node_11(t)*x_node_2(t)