Generalized mass action and S-system power-law models

A biochemical system of power-law rates, in the S-system or the generalized mass action form.

• Name: biochemical_systems_theory
• Field: biology/systems_biology
• Formalism: ordinary differential equations
• References: Raman 2021; Savageau 1976; Voit 2000

Factory

@register_model(
    field='biology/systems_biology',
    tags=('biochemical systems theory', 'power law', 'S-system', 'generalized mass action', 'ode'),
    references=(RAMAN, 'Savageau 1976', 'Voit 2000'),
    defaults={'until': 20.0, 'dt': 0.05},
)
def biochemical_systems_theory(
    x: Sequence[float] = (0.5, 1.0, 4.0),
    *,
    form: str = 's-system',
    alpha: Sequence[float] = PATHWAY_ALPHA,
    beta: Sequence[float] = PATHWAY_BETA,
    g: Sequence[Sequence[float]] = PATHWAY_G,
    h: Sequence[Sequence[float]] = PATHWAY_H,
    stoichiometry: Sequence[Sequence[float]] = PATHWAY_STOICHIOMETRY,
    gamma: Sequence[float] = PATHWAY_GAMMA,
    orders: Sequence[Sequence[float]] = PATHWAY_ORDERS,
) -> System:
    """Return a biochemical system of power-law rates, in the S-system or the generalized mass action form.

    The S-system form gives every species one production and one degradation,
    `dx_i/dt = alpha_i prod_j x_j^g_ij - beta_i prod_j x_j^h_ij`. The generalized mass action form sums power-law
    fluxes, `dx_i/dt = sum_j s_ij gamma_j prod_k x_k^f_jk`. The book states only the two forms. The default is a
    choice, an unbranched pathway `X1 -> X2 -> X3` whose end product inhibits the supply of `X1`, written in both forms
    with the same rates, `dX1/dt = 2 X3^-0.5 - X1^0.5`, `dX2/dt = X1^0.5 - X2^0.5`, and `dX3/dt = X2^0.5 - X3^0.5`.
    The steady state of an S-system solves the linear system `(g - h) log x = log(beta / alpha)`, here `x = (2, 2, 2)`.

    Args:
        x:
            The initial concentrations of the species.

        form:
            `'s-system'` or `'gma'`.

        alpha:
            The rate constants of the productions of the S-system.

        beta:
            The rate constants of the degradations of the S-system.

        g:
            The kinetic orders of the productions of the S-system, one row per species.

        h:
            The kinetic orders of the degradations of the S-system, one row per species.

        stoichiometry:
            The stoichiometric coefficients of the generalized mass action form, one row per species.

        gamma:
            The rate constants of the fluxes of the generalized mass action form.

        orders:
            The kinetic orders of the fluxes of the generalized mass action form, one row per flux.

    Returns:
        system:
            The model, with the concentrations `X1.x`, `X2.x`, and so on.

    Raises:
        ValidationError:
            The form is unknown.

    Examples:
        >>> from skcomplex.models.biology.systems_biology import biochemical_systems_theory
        >>> biochemical_systems_theory(form='gma').interactions[0].gamma
        (2.0, 1.0, 1.0, 1.0)
    """
    species = [SpeciesElement(f'X{index}', x=float(value)) for index, value in enumerate(x, start=1)]
    paths = [item.label for item in species]
    if form == 's-system':
        rates: BaseInteraction = SSystemInteraction(
            'rates',
            species=paths,
            alpha=[float(value) for value in alpha],
            beta=[float(value) for value in beta],
            g=tuple(tuple(row) for row in g),
            h=tuple(tuple(row) for row in h),
        )
    elif form == 'gma':
        rates = GeneralizedMassActionInteraction(
            'rates',
            species=paths,
            stoichiometry=tuple(tuple(row) for row in stoichiometry),
            gamma=tuple(gamma),
            orders=tuple(tuple(row) for row in orders),
        )
    else:
        error_msg = f'Unknown form {form!r}. Use s-system or gma.'
        raise ValidationError(error_msg)
    return System(elements=species, interactions=[rates])

Components

class SpeciesElement(Element):
    """Molecular species, an mRNA, an enzyme, a messenger, or a metabolite, with a nonnegative concentration `x`.

    A species has no laws of its own, and reactions synthesize, convert, and degrade it.

    Args:
        label:
            The label of the species.

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

        x:
            The initial concentration.

    Examples:
        >>> from skcomplex.models.biology.systems_biology import SpeciesElement
        >>> SpeciesElement('camp', x=2.0).x
        2.0
    """

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

class SSystemInteraction(GeneralizedMassActionInteraction):
    """Rates of a biochemical system in the S-system form of biochemical systems theory.

    Every species has one aggregated production and one aggregated degradation, both power laws:
    `dx_i/dt = alpha_i prod_j x_j^g_ij - beta_i prod_j x_j^h_ij`. It is the generalized mass action form with two
    fluxes per species.

    Args:
        label:
            The label of the interaction.

        species:
            The paths of the species.

        alpha:
            The rate constant of the production of every species.

        beta:
            The rate constant of the degradation of every species.

        g:
            The kinetic orders of the productions, one row per species.

        h:
            The kinetic orders of the degradations, one row per species.

    Raises:
        ValidationError:
            The kinetic orders do not have one row and one column per species.

    Examples:
        >>> from skcomplex.models.biology.systems_biology import SSystemInteraction
        >>> SSystemInteraction('rates', species=['X'], alpha=[1.0], beta=[1.0], g=((0.0,),), h=((1.0,),)).h
        ((1.0,),)
    """

    def __init__(
        self,
        label: str,
        *,
        species: Annotated[list[str], RoleSpec(variadic=True)],
        alpha: Annotated[Sequence[float], ParameterSpec(per='species', domain=RealSet.closed_open(0, oo))] = (),
        beta: Annotated[Sequence[float], ParameterSpec(per='species', domain=RealSet.closed_open(0, oo))] = (),
        g: Sequence[Sequence[float]] = ((0.0,),),
        h: Sequence[Sequence[float]] = ((0.0,),),
    ) -> None:
        super().__init__(label, species=species)
        self.alpha = alpha
        self.beta = beta
        self.g = g
        self.h = h

    def _check_matrices(self) -> None:
        """Raise an error for kinetic orders without one row and one column per species."""
        self._check_shape('g', len(self.species), len(self.species))
        self._check_shape('h', len(self.species), len(self.species))

    def _get_fluxes(self, s: SimpleNamespace, values: Sequence[Expression]) -> list[Expression]:
        """Return the production and the degradation flux of every species, in turn."""
        return [
            flux
            for alpha, beta, g, h in zip(s.alpha, s.beta, self.g, self.h, strict=True)
            for flux in (alpha * self._get_power_product(values, g), beta * self._get_power_product(values, h))
        ]

    def _get_stoichiometry(self) -> Sequence[Sequence[float]]:
        """Return the coefficient 1 of every species in its production and -1 in its degradation."""
        size = len(self.species)
        return [
            [1.0 if flux == 2 * row else -1.0 if flux == 2 * row + 1 else 0.0 for flux in range(2 * size)]
            for row in range(size)
        ]

Equations

X1.x: diff(x_X1(t), t) == -1.0*beta_rates_0*sqrt(x_X1(t)) + 1.0*alpha_rates_0/sqrt(x_X3(t))
X2.x: diff(x_X2(t), t) == 1.0*alpha_rates_1*sqrt(x_X1(t)) - 1.0*beta_rates_1*sqrt(x_X2(t))
X3.x: diff(x_X3(t), t) == 1.0*alpha_rates_2*sqrt(x_X2(t)) - 1.0*beta_rates_2*sqrt(x_X3(t))