API Reference

This page contains the auto-generated documentation for the mantis-delta source code.

mantis.CRNetwork

Unified interface for Chemical Reaction Network Theory (CRNT) analysis.

Construct with CRNetwork.from_string(reaction_strings, rates=...). Structural properties (deficiency, conservation laws, weak reversibility) are computed lazily and cached; they do not require rate constants. Symbolic and numerical methods require rate constants to be supplied.

Parameters

reactions : list[Reaction] Directed Reaction objects. Each reversible string A <-> B produces two. rates : dict[str, float], optional Rate constants keyed by reaction strings, e.g. {"A -> B": 1.0}. Keys are automatically normalized to canonical form (species sorted alphabetically); use rate_keys() to inspect the expected form.

Source code in mantis/network.py

class CRNetwork:
    """
    Unified interface for Chemical Reaction Network Theory (CRNT) analysis.

    Construct with `CRNetwork.from_string(reaction_strings, rates=...)`.
    Structural properties (deficiency, conservation laws, weak reversibility) are
    computed lazily and cached; they do not require rate constants.  Symbolic and
    numerical methods require rate constants to be supplied.

    Parameters
    ----------
    reactions : list[Reaction]
        Directed Reaction objects.  Each reversible string `A <-> B` produces two.
    rates : dict[str, float], optional
        Rate constants keyed by reaction strings, e.g. ``{"A -> B": 1.0}``.
        Keys are automatically normalized to canonical form (species sorted
        alphabetically); use ``rate_keys()`` to inspect the expected form.
    """

    def __init__(
        self,
        reactions: list[Reaction],
        rates: dict[str, float] | None = None,
        chemostatted: dict[str, float] | None = None,
    ) -> None:
        self._reactions = reactions
        self._rates: dict[str, float] = {}
        if rates:
            # Normalize user-supplied keys to canonical form
            for key, val in rates.items():
                try:
                    norm = normalize_rate_key(key)
                except ValueError:
                    norm = key
                self._rates[norm] = val
        self._chemostatted: dict[str, float] = dict(chemostatted) if chemostatted else {}

    @classmethod
    def from_string(
        cls,
        reaction_strings: list[str],
        rates: dict[str, float] | None = None,
        chemostatted: dict[str, float] | None = None,
    ) -> CRNetwork:
        """
        Construct from human-readable reaction strings.

        Parameters
        ----------
        reaction_strings : list[str]
            Each string is one reaction, e.g. ``"A + 2B <-> C"`` or ``"ES -> E + P"``.
        rates : dict[str, float], optional
            Maps reaction strings to rate constants.  Order of species within each
            side does not matter; keys are normalized before lookup.
        chemostatted : dict[str, float], optional
            Species held at fixed concentrations by an external reservoir.
            They are excluded from the ODE system and stoichiometry matrix rows,
            but their concentrations are folded into flux expressions.
        """
        reactions = parse_reactions(reaction_strings)
        return cls(reactions, rates, chemostatted)

    # ── Structural properties (no rates needed) ──────────────────────────────

    @cached_property
    def species(self) -> list[str]:
        return build_species_list(self._reactions, set(self._chemostatted.keys()))

    @cached_property
    def complexes(self) -> list[Complex]:
        return list(self._crnt_graph.nodes())

    @cached_property
    def stoichiometry_matrix(self) -> np.ndarray:
        return build_stoichiometry_matrix(self._reactions, self.species)

    @cached_property
    def n_species(self) -> int:
        return len(self.species)

    @cached_property
    def n_complexes(self) -> int:
        return len(self.complexes)

    @cached_property
    def n_reactions(self) -> int:
        return len(self._reactions)

    @cached_property
    def _crnt_graph(self) -> nx.DiGraph:
        return build_complex_graph(self._reactions, set(self._chemostatted.keys()))

    @cached_property
    def n_linkage_classes(self) -> int:
        return len(get_linkage_classes(self._crnt_graph))

    @cached_property
    def deficiency(self) -> int:
        return self._crnt_result.deficiency

    @cached_property
    def is_weakly_reversible(self) -> bool:
        return self._crnt_result.is_weakly_reversible

    @cached_property
    def conservation_laws(self) -> list[sympy.Expr]:
        return conservation_laws_sympy(self.stoichiometry_matrix, self.species)

    @cached_property
    def _crnt_result(self) -> CRNTResult:
        return crnt_analysis(self._reactions, self.stoichiometry_matrix, self.species, set(self._chemostatted.keys()))

    @property
    def chemostatted(self) -> dict[str, float]:
        """Chemostatted species and their fixed concentrations."""
        return dict(self._chemostatted)

    def rate_keys(self) -> list[str]:
        """Return canonical rate key strings for all reactions (for debugging)."""
        return [r.rate_key for r in self._reactions]

    def crnt_summary(self) -> str:
        """Return a human-readable CRNT analysis summary."""
        r = self._crnt_result
        lines = list(r.summary_lines)
        if self._chemostatted:
            chem_str = ", ".join(f"{s}={v}" for s, v in sorted(self._chemostatted.items()))
            lines.insert(0, f"Chemostatted species: {chem_str}")
            lines.insert(1, "")
        lines.append("")
        laws = self.conservation_laws
        lines.append(f"Conservation laws ({len(laws)}):")
        for i, law in enumerate(laws, 1):
            lines.append(f"  [{i}] {law} = const")
        return "\n".join(lines)

    # ── Symbolic ──────────────────────────────────────────────────────────────

    def odes(self, numeric_rates: bool = True) -> dict[str, sympy.Expr]:
        """
        Return mass-action ODEs as SymPy expressions keyed by species name.

        Parameters
        ----------
        numeric_rates : bool
            If True (default) and rates were supplied, substitute numerical values
            so that the returned expressions contain floats rather than ``k_i`` symbols.
        """
        from .symbolic import (
            make_species_symbols,
            make_rate_symbols,
            build_odes,
            substitute_rates,
        )
        sp_syms = make_species_symbols(self.species)
        rate_syms, key_to_sym = make_rate_symbols(self._reactions)
        odes_sym = build_odes(
            self._reactions, self.species, sp_syms, rate_syms, key_to_sym,
            chemostatted_values=self._chemostatted or None,
        )
        if numeric_rates and self._rates:
            odes_sym = {
                sp: substitute_rates(expr, key_to_sym, self._rates)
                for sp, expr in odes_sym.items()
            }
        return odes_sym

    def jacobian(self) -> sympy.Matrix:
        """Return the symbolic Jacobian matrix (n_species × n_species), cached after first call."""
        if "_jacobian_cache" not in self.__dict__:
            from .symbolic import make_species_symbols, build_jacobian
            sp_syms = make_species_symbols(self.species)
            J = build_jacobian(self.odes(numeric_rates=False), self.species, sp_syms)
            self.__dict__["_jacobian_cache"] = J
        return self.__dict__["_jacobian_cache"]

    def jacobian_at(self, ss: dict[str, float]) -> sympy.Matrix:
        from .symbolic import make_species_symbols, substitute_rates, substitute_steady_state
        J_sym = self.jacobian()
        sp_syms = make_species_symbols(self.species)
        from .symbolic import make_rate_symbols
        rate_syms, key_to_sym = make_rate_symbols(self._reactions)
        J_rates = substitute_rates(J_sym, key_to_sym, self._rates)
        return substitute_steady_state(J_rates, sp_syms, ss)

    # ── Numerical ─────────────────────────────────────────────────────────────

    def steady_states(
        self,
        initial_conditions: dict[str, float],
        n_attempts: int = 50,
        seed: int | None = None,
        t_end: float = 1e4,
    ) -> list:
        """
        Find steady states from the given initial conditions.

        Uses ODE integration (Radau) as the primary strategy, which exactly
        preserves conservation laws, then falls back to multi-start least_squares
        for additional steady states. Duplicate solutions (relative L2 distance
        < 10%) and non-physical solutions (any concentration < −tol) are discarded.
        Results are sorted by residual ascending and filtered by relative residual
        against the best solution found to remove spurious stuck-at-IC artifacts.

        Parameters
        ----------
        initial_conditions : dict[str, float]
            Initial concentrations; missing species default to 0.
        n_attempts : int
            Total number of solver attempts (including the integration from IC).
        seed : int, optional
            RNG seed for reproducibility.
        t_end : float
            ODE integration horizon in seconds (default 1e4).  Increase for
            systems with slow reactions whose timescale τ = 1/(k·[X]) >> 1e4 s.

        Returns
        -------
        list[SteadyState]
            Each entry has ``.concentrations``, ``.eigenvalues``, ``.is_stable``,
            ``.is_oscillatory``, and ``.residual``.
        """
        from .analysis import find_steady_states
        return find_steady_states(
            self._reactions,
            self.species,
            self._rates,
            initial_conditions,
            n_attempts=n_attempts,
            seed=seed,
            chemostatted_values=self._chemostatted or None,
            t_end=t_end,
        )

    def all_steady_states(
        self,
        initial_conditions: dict[str, float],
        backend: str = "auto",
        scale: float | None = None,
    ) -> list:
        """
        Exhaustively enumerate *all* positive steady states (with completeness).

        Unlike :meth:`steady_states` (ODE integration + multi-start least-squares,
        which can silently miss unstable fixed points), this solves the mass-action
        steady-state polynomial system for *every* complex root and returns the real
        non-negative ones — so it certifies the full steady-state set rather than
        sampling it.  The classic payoff: on the Schlögl model it returns all three
        equilibria, including the unstable middle branch that integration never
        reaches.

        Parameters
        ----------
        initial_conditions : dict[str, float]
            Used only to fix the conservation-law totals (the stoichiometric
            compatibility class) for closed networks.
        backend : {'auto', 'sympy', 'phcpy'}
            ``'sympy'`` is the pure-Python Gröbner-basis enumerator (default via
            ``'auto'``).  ``'phcpy'`` uses optional polynomial homotopy continuation
            (PHCpack) for larger / stiffer systems — ``pip install mantis-delta[homotopy]``.
        scale : float, optional
            Characteristic concentration used to non-dimensionalise the variables.
            Inferred from the conservation totals / initial conditions when omitted;
            override only if the solver struggles on a badly-scaled network.

        Returns
        -------
        list[SteadyState]
            One entry per distinct positive steady state, sorted by residual, each
            with ``.eigenvalues``, ``.is_stable``, ``.is_oscillatory``.
        """
        from .steady_states import all_steady_states
        return all_steady_states(
            self._reactions,
            self.species,
            self._rates,
            initial_conditions,
            chemostatted_values=self._chemostatted or None,
            backend=backend,
            scale=scale,
        )

    def bifurcation(
        self,
        parameter: str,
        range: tuple[float, float],
        n_points: int = 100,
        initial_conditions: dict[str, float] | None = None,
        plot: bool = False,
        t_end: float = 1e4,
    ) -> object:
        """
        Scan one rate constant over a log-spaced range and collect steady states.

        Parameters
        ----------
        parameter : str
            Rate key to vary, e.g. ``"miR21 + H1 -> miR21_H1"``.
        range : tuple[float, float]
            (min, max) values for the parameter (log scale).
        n_points : int
            Number of points in the scan.
        initial_conditions : dict[str, float], optional
            Initial concentrations for each scan point (defaults to all-zero).
        plot : bool
            If True, display a bifurcation diagram for the first species.
        t_end : float
            ODE integration horizon in seconds (default 1e4).  Increase for
            systems with slow reactions whose timescale τ = 1/(k·[X]) >> 1e4 s.

        Returns
        -------
        BifurcationResult
        """
        from .analysis import scan_bifurcation
        result = scan_bifurcation(
            self._reactions,
            self.species,
            self._rates,
            parameter=parameter,
            param_range=range,
            n_points=n_points,
            initial_conditions=initial_conditions or {},
            n_attempts=20,
            chemostatted_values=self._chemostatted or None,
            t_end=t_end,
        )
        if plot:
            from .plot import plot_bifurcation
            plot_bifurcation(result, species=self.species[0])
        return result

    def continuation(
        self,
        parameter: str,
        range: tuple[float, float],
        initial_conditions: dict[str, float],
        ds: float = 0.05,
        max_steps: int = 2000,
        initial_state: dict[str, float] | None = None,
    ) -> object:
        """
        Trace a steady-state branch by **pseudo-arclength continuation**.

        Unlike :meth:`bifurcation` (an independent log-scan that loses the branch
        at folds and never reports them), this parameterises the branch by
        arclength along the solution curve, so saddle-node folds are traversed
        smoothly and detected.  The classic payoff: on the Schlögl model it traces
        the full S-curve and pins both folds that bound the bistable region.

        Parameters
        ----------
        parameter : str
            Rate key to continue, e.g. ``"A + 2X -> 3X"``.
        range : tuple[float, float]
            ``(λ_min, λ_max)`` box the branch is traced within; to capture an
            S-curve the box must bracket both folds.
        initial_conditions : dict[str, float]
            Fixes the conservation totals (the compatibility class) and seeds the
            branch (a steady state is solved at ``λ_min`` unless ``initial_state``
            is given).
        ds : float
            Arclength step (adapted automatically on Newton failure / success).
        max_steps : int
            Cap on continuation steps (guards against closed isolas).
        initial_state : dict[str, float], optional
            Explicit starting steady state — use to select which branch to follow.

        Returns
        -------
        ContinuationResult
            With ``.parameter_values`` (arclength order — non-monotonic across
            folds), ``.branch``, ``.stable``, ``.bifurcations`` (``.folds()`` /
            ``.hopfs()``), and ``.species_branch(species)``.  ``str()`` summarises.
        """
        from .continuation import pseudo_arclength_continuation
        return pseudo_arclength_continuation(
            self._reactions,
            self.species,
            self._rates,
            parameter=parameter,
            param_range=range,
            initial_conditions=initial_conditions,
            chemostatted_values=self._chemostatted or None,
            ds=ds,
            max_steps=max_steps,
            initial_state=initial_state,
        )

    def simulate(
        self,
        initial_conditions: dict[str, float],
        t_span: tuple[float, float],
        t_eval=None,
        rtol: float = 1e-8,
        atol: float = 1e-12,
    ):
        """
        Integrate the ODE system forward in time and return the full trajectory.

        Parameters
        ----------
        initial_conditions : dict[str, float]
            Initial concentrations; missing species default to 0.
        t_span : (t0, tf)
            Start and end times in seconds.
        t_eval : array-like, optional
            Times at which to store the solution.  Defaults to 200 log-spaced
            points across t_span.
        rtol, atol : float
            Solver tolerances (passed to scipy Radau).

        Returns
        -------
        SimulationResult
            ``.times``, ``.concentrations`` (dict species → 1-D array),
            ``.success``.  Use ``.final()`` for the last time-point dict
            or ``.at(t)`` for a specific time.
        """
        from .analysis import simulate_ode
        return simulate_ode(
            self._reactions,
            self.species,
            self._rates,
            initial_conditions,
            t_span=t_span,
            t_eval=t_eval,
            chemostatted_values=self._chemostatted or None,
            rtol=rtol,
            atol=atol,
        )

    def stochastic_simulate(
        self,
        initial_conditions: dict[str, float],
        t_span: tuple[float, float],
        volume_L: float,
        *,
        initial_as: str = "concentration",
        max_events: int = 1_000_000,
        seed: int | None = None,
    ):
        """
        Single-trajectory Gillespie SSA realization.

        Use when molecule counts are low (∼ ≤ 10³) and the deterministic ODE
        gives the wrong answer — e.g., a CHA cascade at single-cell
        concentrations or stochastic switching in a bistable circuit.

        Parameters
        ----------
        initial_conditions : species → initial count or concentration.
        t_span             : (t0, tf) in seconds.
        volume_L           : reaction volume in liters (e.g. 1e-4 = 100 µL).
        initial_as         : 'concentration' (default) or 'count'.
        max_events         : safety cap on reaction firings.
        seed               : RNG seed for reproducibility.

        Returns
        -------
        StochasticResult with ``.times``, ``.counts``, ``.concentrations``,
        ``.at(t)``, ``.final()``.
        """
        from .analysis import gillespie_simulate
        return gillespie_simulate(
            self._reactions,
            self.species,
            self._rates,
            initial_conditions,
            t_span=t_span,
            volume_L=volume_L,
            initial_as=initial_as,
            max_events=max_events,
            seed=seed,
            chemostatted_values=self._chemostatted or None,
        )

    def stationary_distribution(
        self,
        initial_conditions: dict[str, float],
        volume_L: float,
        *,
        initial_as: str = "concentration",
        max_states: int = 200_000,
    ):
        """
        Exact stochastic stationary distribution (Anderson–Craciun–Kurtz 2010).

        For a complex-balanced network the chemical master equation has a closed-form
        product-of-Poissons stationary distribution with means ``λ_i = N_A·V·c_i*``
        (``c*`` the complex-balanced equilibrium) — written down *from structure
        alone, with no simulation*. With conservation laws it is that product
        conditioned on the conserved totals (e.g. a Binomial for a single moiety).

        Parameters
        ----------
        initial_conditions : dict[str, float]
            Fixes the compatibility class; interpreted per ``initial_as``.
        volume_L : float
            Reaction volume in litres (sets the absolute molecule counts).
        initial_as : 'concentration' or 'count'
        max_states : int
            Cap on enumerated states for the conditioned (closed-network) case.

        Returns
        -------
        StationaryDistribution
            With ``.poisson_means()``, ``.probability(state)``, ``.expected_counts()``,
            ``.marginal(species)`` and ``.sample(...)``.

        Raises
        ------
        ValueError
            If the network is not complex-balanced (the closed form does not apply).
        """
        from .stochastic_stationary import stationary_distribution
        return stationary_distribution(
            self._reactions, self.species, self._rates,
            self.deficiency, self.is_weakly_reversible,
            initial_conditions, volume_L,
            initial_as=initial_as,
            chemostatted_values=self._chemostatted or None,
            max_states=max_states,
        )

    def tau_leap_simulate(
        self,
        initial_conditions: dict[str, float],
        t_span: tuple[float, float],
        volume_L: float,
        *,
        initial_as: str = "concentration",
        tau: float | None = None,
        epsilon: float = 0.03,
        n_record: int = 200,
        seed: int | None = None,
    ):
        """
        τ-leap stochastic simulation (approximate Gillespie).

        Same interface as :meth:`stochastic_simulate` but fires reactions in
        Poisson-distributed bursts over each leap rather than one at a time.
        Roughly N× faster than direct SSA when populations are large (N is
        the mean number of firings per leap), at the cost of asymptotic
        accuracy as populations shrink.

        Parameters
        ----------
        tau : optional fixed leap size (s); if None, adaptive τ per Cao 2006.
        epsilon : adaptive-τ tolerance (max fractional propensity change per
                  leap).  Smaller = more accurate, slower.
        n_record : number of evenly-spaced time points to record.

        See :func:`mantis.analysis.tau_leap_simulate` for full details.
        """
        from .analysis import tau_leap_simulate
        return tau_leap_simulate(
            self._reactions,
            self.species,
            self._rates,
            initial_conditions,
            t_span=t_span,
            volume_L=volume_L,
            initial_as=initial_as,
            tau=tau,
            epsilon=epsilon,
            n_record=n_record,
            seed=seed,
            chemostatted_values=self._chemostatted or None,
        )

    def fsp(
        self,
        initial_conditions: dict[str, float],
        volume_L: float,
        t: float | None = None,
        *,
        initial_as: str = "concentration",
        max_states: int = 100_000,
    ):
        """
        Solve the chemical master equation by Finite State Projection (Munsky–Khammash 2006).

        Truncates the count-state space to the reachable set, builds the CME generator,
        and solves it directly — accurate where the SSA is noisy (low molecule numbers),
        with a rigorous truncation-error bound.

        Parameters
        ----------
        initial_conditions : dict[str, float]
            Starting state; interpreted per ``initial_as``.
        volume_L : float
            Reaction volume in litres.
        t : float, optional
            If given, returns the transient distribution ``p(t)``; if ``None``, returns
            the stationary distribution (valid for closed, fully enumerated networks).
        initial_as : 'concentration' or 'count'
        max_states : int
            Cap on the number of enumerated states.

        Returns
        -------
        FSPResult
            With ``.probabilities``, ``.states``, ``.marginal(species)``,
            ``.expected_counts()``, and ``.truncation_error``.
        """
        from .fsp import fsp_solve
        return fsp_solve(
            self._reactions, self.species, self._rates,
            initial_conditions, volume_L, t=t, initial_as=initial_as,
            chemostatted_values=self._chemostatted or None,
            max_states=max_states,
        )

    def draw(self, ax=None, layout: str = "spring") -> matplotlib.axes.Axes:
        from .plot import draw_reaction_graph
        return draw_reaction_graph(self._crnt_graph, ax=ax, layout=layout)

    # ── Modern structural analysis ──────────────────────────────────────────────

    def detect_acr(
        self, initial_conditions: dict[str, float] | None = None
    ):
        """
        Detect Absolute Concentration Robustness (Shinar–Feinberg, *Science* 2010).

        A deficiency-one network with two non-terminal complexes differing in a single
        species *S* has ACR in *S*: its steady-state concentration is independent of the
        initial totals — the robustness property at the heart of biosensor design.

        Parameters
        ----------
        initial_conditions : dict[str, float], optional
            When supplied (with rate constants set), the robust steady-state value of
            each ACR species is computed at a positive steady state.

        Returns
        -------
        ACRResult
            With ``.has_acr``, ``.species``, ``.acr_values``, and a readable summary
            via ``str()``.
        """
        from .acr import detect_acr
        return detect_acr(
            self._reactions, self.species, self._crnt_graph, self.deficiency,
            rate_values=self._rates or None,
            initial_conditions=initial_conditions,
            chemostatted_values=self._chemostatted or None,
        )

    def is_injective(self):
        """
        Test mass-action injectivity (Craciun–Feinberg, *SIAM J. Appl. Math.* 2005/06).

        An injective network has **at most one** positive steady state in each
        stoichiometric compatibility class — ruling out multistationarity. Unlike the
        Deficiency One Theorem this works for higher-deficiency networks too. The test
        is the sign-definiteness of the steady-state Jacobian determinant over all
        positive concentrations and rate constants (a sufficient condition).

        Returns
        -------
        InjectivityResult
            With ``.injective`` (True ⇒ certified injective; False ⇒ inconclusive),
            ``.determinant``, and a readable summary via ``str()``.
        """
        from .injectivity import test_injectivity
        return test_injectivity(
            self._reactions, self.species,
            chemostatted_values=self._chemostatted or None,
        )

    def multistationarity_region(self):
        """
        Compute the multistationarity parameter region (Conradi–Feliu–Mincheva–Wiuf,
        *PLoS Comput. Biol.* 2017).

        Where :meth:`is_injective` answers *whether* multistationarity is possible,
        this maps *where in parameter space* it occurs.  It forms the critical
        function ``φ = (−1)^s·det J`` (``s = rank(N)``) and reads it as a polynomial
        in the steady-state concentrations: if every coefficient is positive for all
        positive rates the network is monostationary for all parameters; otherwise
        each coefficient that can go negative is a face of the multistationarity
        region.  When rate constants are set, also reports whether *those* rates lie
        inside the region.

        Returns
        -------
        MultistationarityResult
            With ``.monostationary``, ``.multistationary_possible``,
            ``.critical_function``, ``.region_conditions``,
            ``.multistationary_at_rates``, and a readable summary via ``str()``.
        """
        from .multistationarity import multistationarity_region
        return multistationarity_region(
            self._reactions, self.species,
            rate_values=self._rates or None,
            chemostatted_values=self._chemostatted or None,
        )

    def is_multistationary(
        self, initial_conditions: dict[str, float], backend: str = "auto"
    ) -> bool:
        """
        Exactly decide multistationarity in one compatibility class by enumeration.

        Whereas :meth:`multistationarity_region` gives the rate-free *structural*
        verdict (and the symbolic region), this is the concrete numeric answer for
        the supplied rate constants and totals: it enumerates *all* positive steady
        states (:meth:`all_steady_states`) in the class fixed by ``initial_conditions``
        and returns ``True`` iff there are at least two.  Definitive for any network —
        no parametrisation or sign condition required.

        Parameters
        ----------
        initial_conditions : dict[str, float]
            Fixes the stoichiometric compatibility class (the conservation totals).
        backend : {'auto', 'sympy', 'phcpy'}
            Steady-state enumeration backend (see :meth:`all_steady_states`).

        Returns
        -------
        bool
            ``True`` ⇔ ≥ 2 positive steady states in this class.
        """
        states = self.all_steady_states(initial_conditions, backend=backend)
        return len(states) >= 2

    # ── Global stability (Horn–Jackson Lyapunov) ────────────────────────────────

    def is_complex_balanced(
        self, initial_conditions: dict[str, float] | None = None
    ) -> bool:
        """
        Whether the network is complex-balanced.

        Every weakly reversible deficiency-zero network is complex-balanced for
        *all* rate constants (Horn–Jackson). For weakly reversible higher-deficiency
        networks the property is rate-dependent and is verified numerically at the
        positive equilibrium — pass ``initial_conditions`` to enable that check.
        """
        from .stability import is_complex_balanced
        return is_complex_balanced(
            self._reactions, self.species, self._rates,
            self.deficiency, self.is_weakly_reversible,
            initial_conditions=initial_conditions,
            chemostatted_values=self._chemostatted or None,
        )

    def complex_balanced_equilibrium(
        self, initial_conditions: dict[str, float]
    ) -> dict[str, float]:
        """
        Return the **Birch point** — the unique positive complex-balanced equilibrium
        in the stoichiometric compatibility class fixed by ``initial_conditions``.

        Defined only for complex-balanced networks (every weakly reversible
        deficiency-zero network qualifies, for all rates). This equilibrium sets the
        means of the exact stochastic stationary distribution (see
        :meth:`stationary_distribution`). Raises ``ValueError`` otherwise.
        """
        from .stability import complex_balanced_equilibrium
        return complex_balanced_equilibrium(
            self._reactions, self.species, self._rates,
            self.deficiency, self.is_weakly_reversible, initial_conditions,
            chemostatted_values=self._chemostatted or None,
        )

    def certify_global_stability(
        self,
        initial_conditions: dict[str, float],
        n_samples: int = 200,
        seed: int = 0,
    ):
        """
        Certify *global* asymptotic stability via the Horn–Jackson Lyapunov function.

        For complex-balanced networks the pseudo-Helmholtz function
        ``V(c) = Σ_i [c_i(ln(c_i/c_i*) − 1) + c_i*]`` is a strict Lyapunov function,
        so the positive equilibrium ``c*`` is globally asymptotically stable within
        its stoichiometric compatibility class — a strictly stronger statement than
        the local eigenvalue test.

        Parameters
        ----------
        initial_conditions : dict[str, float]
            Fixes the compatibility class (conservation totals) and the equilibrium.
        n_samples : int
            Number of random points in the class at which ``dV/dt ≤ 0`` is verified.
        seed : int
            RNG seed for the sampling.

        Returns
        -------
        StabilityCertificate
            With ``.globally_stable``, ``.equilibrium``, ``.lyapunov_function`` (V),
            and ``.lyapunov_derivative`` (dV/dt).  ``str()`` gives a readable summary.
        """
        from .stability import certify_global_stability
        return certify_global_stability(
            self._reactions, self.species, self._rates,
            self.deficiency, self.is_weakly_reversible,
            self.odes(numeric_rates=True), initial_conditions,
            chemostatted_values=self._chemostatted or None,
            n_samples=n_samples, seed=seed,
        )

    # ── Interoperability (SBML) ─────────────────────────────────────────────────

    def to_sbml(self, path: str | None = None) -> str:
        """
        Export this network to SBML (Level 3 Version 2).

        Requires the optional ``python-libsbml`` package
        (``pip install mantis-delta[sbml]``).

        Parameters
        ----------
        path : str, optional
            If given, the SBML document is written to this file path.  The XML
            string is returned in all cases.

        Returns
        -------
        str
            The SBML document as an XML string.
        """
        from .sbml import network_to_sbml, write_sbml
        if path is not None:
            write_sbml(self, path)
        return network_to_sbml(self)

    @classmethod
    def from_sbml(cls, source: str) -> CRNetwork:
        """
        Construct a network from an SBML model.

        Requires the optional ``python-libsbml`` package
        (``pip install mantis-delta[sbml]``).

        Parameters
        ----------
        source : str
            Either a path to an ``.xml`` SBML file or a raw SBML/XML string.
            Mass-action kinetic laws are parsed to recover rate constants where
            possible; SBML boundary/constant species become chemostatted species.

        Returns
        -------
        CRNetwork
        """
        from .sbml import network_from_sbml
        return network_from_sbml(source)

chemostatted property

Chemostatted species and their fixed concentrations.

all_steady_states(initial_conditions, backend='auto', scale=None)

Exhaustively enumerate all positive steady states (with completeness).

Unlike :meth:steady_states (ODE integration + multi-start least-squares, which can silently miss unstable fixed points), this solves the mass-action steady-state polynomial system for every complex root and returns the real non-negative ones — so it certifies the full steady-state set rather than sampling it. The classic payoff: on the Schlögl model it returns all three equilibria, including the unstable middle branch that integration never reaches.

Parameters

initial_conditions : dict[str, float] Used only to fix the conservation-law totals (the stoichiometric compatibility class) for closed networks. backend : {'auto', 'sympy', 'phcpy'} 'sympy' is the pure-Python Gröbner-basis enumerator (default via 'auto'). 'phcpy' uses optional polynomial homotopy continuation (PHCpack) for larger / stiffer systems — pip install mantis-delta[homotopy]. scale : float, optional Characteristic concentration used to non-dimensionalise the variables. Inferred from the conservation totals / initial conditions when omitted; override only if the solver struggles on a badly-scaled network.

Returns

list[SteadyState] One entry per distinct positive steady state, sorted by residual, each with .eigenvalues, .is_stable, .is_oscillatory.

Source code in mantis/network.py

def all_steady_states(
    self,
    initial_conditions: dict[str, float],
    backend: str = "auto",
    scale: float | None = None,
) -> list:
    """
    Exhaustively enumerate *all* positive steady states (with completeness).

    Unlike :meth:`steady_states` (ODE integration + multi-start least-squares,
    which can silently miss unstable fixed points), this solves the mass-action
    steady-state polynomial system for *every* complex root and returns the real
    non-negative ones — so it certifies the full steady-state set rather than
    sampling it.  The classic payoff: on the Schlögl model it returns all three
    equilibria, including the unstable middle branch that integration never
    reaches.

    Parameters
    ----------
    initial_conditions : dict[str, float]
        Used only to fix the conservation-law totals (the stoichiometric
        compatibility class) for closed networks.
    backend : {'auto', 'sympy', 'phcpy'}
        ``'sympy'`` is the pure-Python Gröbner-basis enumerator (default via
        ``'auto'``).  ``'phcpy'`` uses optional polynomial homotopy continuation
        (PHCpack) for larger / stiffer systems — ``pip install mantis-delta[homotopy]``.
    scale : float, optional
        Characteristic concentration used to non-dimensionalise the variables.
        Inferred from the conservation totals / initial conditions when omitted;
        override only if the solver struggles on a badly-scaled network.

    Returns
    -------
    list[SteadyState]
        One entry per distinct positive steady state, sorted by residual, each
        with ``.eigenvalues``, ``.is_stable``, ``.is_oscillatory``.
    """
    from .steady_states import all_steady_states
    return all_steady_states(
        self._reactions,
        self.species,
        self._rates,
        initial_conditions,
        chemostatted_values=self._chemostatted or None,
        backend=backend,
        scale=scale,
    )

bifurcation(parameter, range, n_points=100, initial_conditions=None, plot=False, t_end=10000.0)

Scan one rate constant over a log-spaced range and collect steady states.

Parameters

parameter : str Rate key to vary, e.g. "miR21 + H1 -> miR21_H1". range : tuple[float, float] (min, max) values for the parameter (log scale). n_points : int Number of points in the scan. initial_conditions : dict[str, float], optional Initial concentrations for each scan point (defaults to all-zero). plot : bool If True, display a bifurcation diagram for the first species. t_end : float ODE integration horizon in seconds (default 1e4). Increase for systems with slow reactions whose timescale τ = 1/(k·[X]) >> 1e4 s.

Returns

BifurcationResult

Source code in mantis/network.py

def bifurcation(
    self,
    parameter: str,
    range: tuple[float, float],
    n_points: int = 100,
    initial_conditions: dict[str, float] | None = None,
    plot: bool = False,
    t_end: float = 1e4,
) -> object:
    """
    Scan one rate constant over a log-spaced range and collect steady states.

    Parameters
    ----------
    parameter : str
        Rate key to vary, e.g. ``"miR21 + H1 -> miR21_H1"``.
    range : tuple[float, float]
        (min, max) values for the parameter (log scale).
    n_points : int
        Number of points in the scan.
    initial_conditions : dict[str, float], optional
        Initial concentrations for each scan point (defaults to all-zero).
    plot : bool
        If True, display a bifurcation diagram for the first species.
    t_end : float
        ODE integration horizon in seconds (default 1e4).  Increase for
        systems with slow reactions whose timescale τ = 1/(k·[X]) >> 1e4 s.

    Returns
    -------
    BifurcationResult
    """
    from .analysis import scan_bifurcation
    result = scan_bifurcation(
        self._reactions,
        self.species,
        self._rates,
        parameter=parameter,
        param_range=range,
        n_points=n_points,
        initial_conditions=initial_conditions or {},
        n_attempts=20,
        chemostatted_values=self._chemostatted or None,
        t_end=t_end,
    )
    if plot:
        from .plot import plot_bifurcation
        plot_bifurcation(result, species=self.species[0])
    return result

certify_global_stability(initial_conditions, n_samples=200, seed=0)

Certify global asymptotic stability via the Horn–Jackson Lyapunov function.

For complex-balanced networks the pseudo-Helmholtz function V(c) = Σ_i [c_i(ln(c_i/c_i*) − 1) + c_i*] is a strict Lyapunov function, so the positive equilibrium c* is globally asymptotically stable within its stoichiometric compatibility class — a strictly stronger statement than the local eigenvalue test.

Parameters

initial_conditions : dict[str, float] Fixes the compatibility class (conservation totals) and the equilibrium. n_samples : int Number of random points in the class at which dV/dt ≤ 0 is verified. seed : int RNG seed for the sampling.

Returns

StabilityCertificate With .globally_stable, .equilibrium, .lyapunov_function (V), and .lyapunov_derivative (dV/dt). str() gives a readable summary.

Source code in mantis/network.py

def certify_global_stability(
    self,
    initial_conditions: dict[str, float],
    n_samples: int = 200,
    seed: int = 0,
):
    """
    Certify *global* asymptotic stability via the Horn–Jackson Lyapunov function.

    For complex-balanced networks the pseudo-Helmholtz function
    ``V(c) = Σ_i [c_i(ln(c_i/c_i*) − 1) + c_i*]`` is a strict Lyapunov function,
    so the positive equilibrium ``c*`` is globally asymptotically stable within
    its stoichiometric compatibility class — a strictly stronger statement than
    the local eigenvalue test.

    Parameters
    ----------
    initial_conditions : dict[str, float]
        Fixes the compatibility class (conservation totals) and the equilibrium.
    n_samples : int
        Number of random points in the class at which ``dV/dt ≤ 0`` is verified.
    seed : int
        RNG seed for the sampling.

    Returns
    -------
    StabilityCertificate
        With ``.globally_stable``, ``.equilibrium``, ``.lyapunov_function`` (V),
        and ``.lyapunov_derivative`` (dV/dt).  ``str()`` gives a readable summary.
    """
    from .stability import certify_global_stability
    return certify_global_stability(
        self._reactions, self.species, self._rates,
        self.deficiency, self.is_weakly_reversible,
        self.odes(numeric_rates=True), initial_conditions,
        chemostatted_values=self._chemostatted or None,
        n_samples=n_samples, seed=seed,
    )

complex_balanced_equilibrium(initial_conditions)

Return the Birch point — the unique positive complex-balanced equilibrium in the stoichiometric compatibility class fixed by initial_conditions.

Defined only for complex-balanced networks (every weakly reversible deficiency-zero network qualifies, for all rates). This equilibrium sets the means of the exact stochastic stationary distribution (see :meth:stationary_distribution). Raises ValueError otherwise.

Source code in mantis/network.py

def complex_balanced_equilibrium(
    self, initial_conditions: dict[str, float]
) -> dict[str, float]:
    """
    Return the **Birch point** — the unique positive complex-balanced equilibrium
    in the stoichiometric compatibility class fixed by ``initial_conditions``.

    Defined only for complex-balanced networks (every weakly reversible
    deficiency-zero network qualifies, for all rates). This equilibrium sets the
    means of the exact stochastic stationary distribution (see
    :meth:`stationary_distribution`). Raises ``ValueError`` otherwise.
    """
    from .stability import complex_balanced_equilibrium
    return complex_balanced_equilibrium(
        self._reactions, self.species, self._rates,
        self.deficiency, self.is_weakly_reversible, initial_conditions,
        chemostatted_values=self._chemostatted or None,
    )

continuation(parameter, range, initial_conditions, ds=0.05, max_steps=2000, initial_state=None)

Trace a steady-state branch by pseudo-arclength continuation.

Unlike :meth:bifurcation (an independent log-scan that loses the branch at folds and never reports them), this parameterises the branch by arclength along the solution curve, so saddle-node folds are traversed smoothly and detected. The classic payoff: on the Schlögl model it traces the full S-curve and pins both folds that bound the bistable region.

Parameters

parameter : str Rate key to continue, e.g. "A + 2X -> 3X". range : tuple[float, float] (λ_min, λ_max) box the branch is traced within; to capture an S-curve the box must bracket both folds. initial_conditions : dict[str, float] Fixes the conservation totals (the compatibility class) and seeds the branch (a steady state is solved at λ_min unless initial_state is given). ds : float Arclength step (adapted automatically on Newton failure / success). max_steps : int Cap on continuation steps (guards against closed isolas). initial_state : dict[str, float], optional Explicit starting steady state — use to select which branch to follow.

Returns

ContinuationResult With .parameter_values (arclength order — non-monotonic across folds), .branch, .stable, .bifurcations (.folds() / .hopfs()), and .species_branch(species). str() summarises.

Source code in mantis/network.py

def continuation(
    self,
    parameter: str,
    range: tuple[float, float],
    initial_conditions: dict[str, float],
    ds: float = 0.05,
    max_steps: int = 2000,
    initial_state: dict[str, float] | None = None,
) -> object:
    """
    Trace a steady-state branch by **pseudo-arclength continuation**.

    Unlike :meth:`bifurcation` (an independent log-scan that loses the branch
    at folds and never reports them), this parameterises the branch by
    arclength along the solution curve, so saddle-node folds are traversed
    smoothly and detected.  The classic payoff: on the Schlögl model it traces
    the full S-curve and pins both folds that bound the bistable region.

    Parameters
    ----------
    parameter : str
        Rate key to continue, e.g. ``"A + 2X -> 3X"``.
    range : tuple[float, float]
        ``(λ_min, λ_max)`` box the branch is traced within; to capture an
        S-curve the box must bracket both folds.
    initial_conditions : dict[str, float]
        Fixes the conservation totals (the compatibility class) and seeds the
        branch (a steady state is solved at ``λ_min`` unless ``initial_state``
        is given).
    ds : float
        Arclength step (adapted automatically on Newton failure / success).
    max_steps : int
        Cap on continuation steps (guards against closed isolas).
    initial_state : dict[str, float], optional
        Explicit starting steady state — use to select which branch to follow.

    Returns
    -------
    ContinuationResult
        With ``.parameter_values`` (arclength order — non-monotonic across
        folds), ``.branch``, ``.stable``, ``.bifurcations`` (``.folds()`` /
        ``.hopfs()``), and ``.species_branch(species)``.  ``str()`` summarises.
    """
    from .continuation import pseudo_arclength_continuation
    return pseudo_arclength_continuation(
        self._reactions,
        self.species,
        self._rates,
        parameter=parameter,
        param_range=range,
        initial_conditions=initial_conditions,
        chemostatted_values=self._chemostatted or None,
        ds=ds,
        max_steps=max_steps,
        initial_state=initial_state,
    )

crnt_summary()

Return a human-readable CRNT analysis summary.

Source code in mantis/network.py

def crnt_summary(self) -> str:
    """Return a human-readable CRNT analysis summary."""
    r = self._crnt_result
    lines = list(r.summary_lines)
    if self._chemostatted:
        chem_str = ", ".join(f"{s}={v}" for s, v in sorted(self._chemostatted.items()))
        lines.insert(0, f"Chemostatted species: {chem_str}")
        lines.insert(1, "")
    lines.append("")
    laws = self.conservation_laws
    lines.append(f"Conservation laws ({len(laws)}):")
    for i, law in enumerate(laws, 1):
        lines.append(f"  [{i}] {law} = const")
    return "\n".join(lines)

detect_acr(initial_conditions=None)

Detect Absolute Concentration Robustness (Shinar–Feinberg, Science 2010).

A deficiency-one network with two non-terminal complexes differing in a single species S has ACR in S: its steady-state concentration is independent of the initial totals — the robustness property at the heart of biosensor design.

Parameters

initial_conditions : dict[str, float], optional When supplied (with rate constants set), the robust steady-state value of each ACR species is computed at a positive steady state.

Returns

ACRResult With .has_acr, .species, .acr_values, and a readable summary via str().

Source code in mantis/network.py

def detect_acr(
    self, initial_conditions: dict[str, float] | None = None
):
    """
    Detect Absolute Concentration Robustness (Shinar–Feinberg, *Science* 2010).

    A deficiency-one network with two non-terminal complexes differing in a single
    species *S* has ACR in *S*: its steady-state concentration is independent of the
    initial totals — the robustness property at the heart of biosensor design.

    Parameters
    ----------
    initial_conditions : dict[str, float], optional
        When supplied (with rate constants set), the robust steady-state value of
        each ACR species is computed at a positive steady state.

    Returns
    -------
    ACRResult
        With ``.has_acr``, ``.species``, ``.acr_values``, and a readable summary
        via ``str()``.
    """
    from .acr import detect_acr
    return detect_acr(
        self._reactions, self.species, self._crnt_graph, self.deficiency,
        rate_values=self._rates or None,
        initial_conditions=initial_conditions,
        chemostatted_values=self._chemostatted or None,
    )

from_sbml(source) classmethod

Construct a network from an SBML model.

Requires the optional python-libsbml package (pip install mantis-delta[sbml]).

Parameters

source : str Either a path to an .xml SBML file or a raw SBML/XML string. Mass-action kinetic laws are parsed to recover rate constants where possible; SBML boundary/constant species become chemostatted species.

Returns

CRNetwork

Source code in mantis/network.py

@classmethod
def from_sbml(cls, source: str) -> CRNetwork:
    """
    Construct a network from an SBML model.

    Requires the optional ``python-libsbml`` package
    (``pip install mantis-delta[sbml]``).

    Parameters
    ----------
    source : str
        Either a path to an ``.xml`` SBML file or a raw SBML/XML string.
        Mass-action kinetic laws are parsed to recover rate constants where
        possible; SBML boundary/constant species become chemostatted species.

    Returns
    -------
    CRNetwork
    """
    from .sbml import network_from_sbml
    return network_from_sbml(source)

from_string(reaction_strings, rates=None, chemostatted=None) classmethod

Construct from human-readable reaction strings.

Parameters

reaction_strings : list[str] Each string is one reaction, e.g. "A + 2B <-> C" or "ES -> E + P". rates : dict[str, float], optional Maps reaction strings to rate constants. Order of species within each side does not matter; keys are normalized before lookup. chemostatted : dict[str, float], optional Species held at fixed concentrations by an external reservoir. They are excluded from the ODE system and stoichiometry matrix rows, but their concentrations are folded into flux expressions.

Source code in mantis/network.py

@classmethod
def from_string(
    cls,
    reaction_strings: list[str],
    rates: dict[str, float] | None = None,
    chemostatted: dict[str, float] | None = None,
) -> CRNetwork:
    """
    Construct from human-readable reaction strings.

    Parameters
    ----------
    reaction_strings : list[str]
        Each string is one reaction, e.g. ``"A + 2B <-> C"`` or ``"ES -> E + P"``.
    rates : dict[str, float], optional
        Maps reaction strings to rate constants.  Order of species within each
        side does not matter; keys are normalized before lookup.
    chemostatted : dict[str, float], optional
        Species held at fixed concentrations by an external reservoir.
        They are excluded from the ODE system and stoichiometry matrix rows,
        but their concentrations are folded into flux expressions.
    """
    reactions = parse_reactions(reaction_strings)
    return cls(reactions, rates, chemostatted)

fsp(initial_conditions, volume_L, t=None, *, initial_as='concentration', max_states=100000)

Solve the chemical master equation by Finite State Projection (Munsky–Khammash 2006).

Truncates the count-state space to the reachable set, builds the CME generator, and solves it directly — accurate where the SSA is noisy (low molecule numbers), with a rigorous truncation-error bound.

Parameters

initial_conditions : dict[str, float] Starting state; interpreted per initial_as. volume_L : float Reaction volume in litres. t : float, optional If given, returns the transient distribution p(t); if None, returns the stationary distribution (valid for closed, fully enumerated networks). initial_as : 'concentration' or 'count' max_states : int Cap on the number of enumerated states.

Returns

FSPResult With .probabilities, .states, .marginal(species), .expected_counts(), and .truncation_error.

Source code in mantis/network.py

def fsp(
    self,
    initial_conditions: dict[str, float],
    volume_L: float,
    t: float | None = None,
    *,
    initial_as: str = "concentration",
    max_states: int = 100_000,
):
    """
    Solve the chemical master equation by Finite State Projection (Munsky–Khammash 2006).

    Truncates the count-state space to the reachable set, builds the CME generator,
    and solves it directly — accurate where the SSA is noisy (low molecule numbers),
    with a rigorous truncation-error bound.

    Parameters
    ----------
    initial_conditions : dict[str, float]
        Starting state; interpreted per ``initial_as``.
    volume_L : float
        Reaction volume in litres.
    t : float, optional
        If given, returns the transient distribution ``p(t)``; if ``None``, returns
        the stationary distribution (valid for closed, fully enumerated networks).
    initial_as : 'concentration' or 'count'
    max_states : int
        Cap on the number of enumerated states.

    Returns
    -------
    FSPResult
        With ``.probabilities``, ``.states``, ``.marginal(species)``,
        ``.expected_counts()``, and ``.truncation_error``.
    """
    from .fsp import fsp_solve
    return fsp_solve(
        self._reactions, self.species, self._rates,
        initial_conditions, volume_L, t=t, initial_as=initial_as,
        chemostatted_values=self._chemostatted or None,
        max_states=max_states,
    )

is_complex_balanced(initial_conditions=None)

Whether the network is complex-balanced.

Every weakly reversible deficiency-zero network is complex-balanced for all rate constants (Horn–Jackson). For weakly reversible higher-deficiency networks the property is rate-dependent and is verified numerically at the positive equilibrium — pass initial_conditions to enable that check.

Source code in mantis/network.py

def is_complex_balanced(
    self, initial_conditions: dict[str, float] | None = None
) -> bool:
    """
    Whether the network is complex-balanced.

    Every weakly reversible deficiency-zero network is complex-balanced for
    *all* rate constants (Horn–Jackson). For weakly reversible higher-deficiency
    networks the property is rate-dependent and is verified numerically at the
    positive equilibrium — pass ``initial_conditions`` to enable that check.
    """
    from .stability import is_complex_balanced
    return is_complex_balanced(
        self._reactions, self.species, self._rates,
        self.deficiency, self.is_weakly_reversible,
        initial_conditions=initial_conditions,
        chemostatted_values=self._chemostatted or None,
    )

is_injective()

Test mass-action injectivity (Craciun–Feinberg, SIAM J. Appl. Math. 2005/06).

An injective network has at most one positive steady state in each stoichiometric compatibility class — ruling out multistationarity. Unlike the Deficiency One Theorem this works for higher-deficiency networks too. The test is the sign-definiteness of the steady-state Jacobian determinant over all positive concentrations and rate constants (a sufficient condition).

Returns

InjectivityResult With .injective (True ⇒ certified injective; False ⇒ inconclusive), .determinant, and a readable summary via str().

Source code in mantis/network.py

def is_injective(self):
    """
    Test mass-action injectivity (Craciun–Feinberg, *SIAM J. Appl. Math.* 2005/06).

    An injective network has **at most one** positive steady state in each
    stoichiometric compatibility class — ruling out multistationarity. Unlike the
    Deficiency One Theorem this works for higher-deficiency networks too. The test
    is the sign-definiteness of the steady-state Jacobian determinant over all
    positive concentrations and rate constants (a sufficient condition).

    Returns
    -------
    InjectivityResult
        With ``.injective`` (True ⇒ certified injective; False ⇒ inconclusive),
        ``.determinant``, and a readable summary via ``str()``.
    """
    from .injectivity import test_injectivity
    return test_injectivity(
        self._reactions, self.species,
        chemostatted_values=self._chemostatted or None,
    )

is_multistationary(initial_conditions, backend='auto')

Exactly decide multistationarity in one compatibility class by enumeration.

Whereas :meth:multistationarity_region gives the rate-free structural verdict (and the symbolic region), this is the concrete numeric answer for the supplied rate constants and totals: it enumerates all positive steady states (:meth:all_steady_states) in the class fixed by initial_conditions and returns True iff there are at least two. Definitive for any network — no parametrisation or sign condition required.

Parameters

initial_conditions : dict[str, float] Fixes the stoichiometric compatibility class (the conservation totals). backend : {'auto', 'sympy', 'phcpy'} Steady-state enumeration backend (see :meth:all_steady_states).

Returns

bool True ⇔ ≥ 2 positive steady states in this class.

Source code in mantis/network.py

def is_multistationary(
    self, initial_conditions: dict[str, float], backend: str = "auto"
) -> bool:
    """
    Exactly decide multistationarity in one compatibility class by enumeration.

    Whereas :meth:`multistationarity_region` gives the rate-free *structural*
    verdict (and the symbolic region), this is the concrete numeric answer for
    the supplied rate constants and totals: it enumerates *all* positive steady
    states (:meth:`all_steady_states`) in the class fixed by ``initial_conditions``
    and returns ``True`` iff there are at least two.  Definitive for any network —
    no parametrisation or sign condition required.

    Parameters
    ----------
    initial_conditions : dict[str, float]
        Fixes the stoichiometric compatibility class (the conservation totals).
    backend : {'auto', 'sympy', 'phcpy'}
        Steady-state enumeration backend (see :meth:`all_steady_states`).

    Returns
    -------
    bool
        ``True`` ⇔ ≥ 2 positive steady states in this class.
    """
    states = self.all_steady_states(initial_conditions, backend=backend)
    return len(states) >= 2

jacobian()

Return the symbolic Jacobian matrix (n_species × n_species), cached after first call.

Source code in mantis/network.py

def jacobian(self) -> sympy.Matrix:
    """Return the symbolic Jacobian matrix (n_species × n_species), cached after first call."""
    if "_jacobian_cache" not in self.__dict__:
        from .symbolic import make_species_symbols, build_jacobian
        sp_syms = make_species_symbols(self.species)
        J = build_jacobian(self.odes(numeric_rates=False), self.species, sp_syms)
        self.__dict__["_jacobian_cache"] = J
    return self.__dict__["_jacobian_cache"]

multistationarity_region()

Compute the multistationarity parameter region (Conradi–Feliu–Mincheva–Wiuf, PLoS Comput. Biol. 2017).

Where :meth:is_injective answers whether multistationarity is possible, this maps where in parameter space it occurs. It forms the critical function φ = (−1)^s·det J (s = rank(N)) and reads it as a polynomial in the steady-state concentrations: if every coefficient is positive for all positive rates the network is monostationary for all parameters; otherwise each coefficient that can go negative is a face of the multistationarity region. When rate constants are set, also reports whether those rates lie inside the region.

Returns

MultistationarityResult With .monostationary, .multistationary_possible, .critical_function, .region_conditions, .multistationary_at_rates, and a readable summary via str().

Source code in mantis/network.py

def multistationarity_region(self):
    """
    Compute the multistationarity parameter region (Conradi–Feliu–Mincheva–Wiuf,
    *PLoS Comput. Biol.* 2017).

    Where :meth:`is_injective` answers *whether* multistationarity is possible,
    this maps *where in parameter space* it occurs.  It forms the critical
    function ``φ = (−1)^s·det J`` (``s = rank(N)``) and reads it as a polynomial
    in the steady-state concentrations: if every coefficient is positive for all
    positive rates the network is monostationary for all parameters; otherwise
    each coefficient that can go negative is a face of the multistationarity
    region.  When rate constants are set, also reports whether *those* rates lie
    inside the region.

    Returns
    -------
    MultistationarityResult
        With ``.monostationary``, ``.multistationary_possible``,
        ``.critical_function``, ``.region_conditions``,
        ``.multistationary_at_rates``, and a readable summary via ``str()``.
    """
    from .multistationarity import multistationarity_region
    return multistationarity_region(
        self._reactions, self.species,
        rate_values=self._rates or None,
        chemostatted_values=self._chemostatted or None,
    )

odes(numeric_rates=True)

Return mass-action ODEs as SymPy expressions keyed by species name.

Parameters

numeric_rates : bool If True (default) and rates were supplied, substitute numerical values so that the returned expressions contain floats rather than k_i symbols.

Source code in mantis/network.py

def odes(self, numeric_rates: bool = True) -> dict[str, sympy.Expr]:
    """
    Return mass-action ODEs as SymPy expressions keyed by species name.

    Parameters
    ----------
    numeric_rates : bool
        If True (default) and rates were supplied, substitute numerical values
        so that the returned expressions contain floats rather than ``k_i`` symbols.
    """
    from .symbolic import (
        make_species_symbols,
        make_rate_symbols,
        build_odes,
        substitute_rates,
    )
    sp_syms = make_species_symbols(self.species)
    rate_syms, key_to_sym = make_rate_symbols(self._reactions)
    odes_sym = build_odes(
        self._reactions, self.species, sp_syms, rate_syms, key_to_sym,
        chemostatted_values=self._chemostatted or None,
    )
    if numeric_rates and self._rates:
        odes_sym = {
            sp: substitute_rates(expr, key_to_sym, self._rates)
            for sp, expr in odes_sym.items()
        }
    return odes_sym

rate_keys()

Return canonical rate key strings for all reactions (for debugging).

Source code in mantis/network.py

def rate_keys(self) -> list[str]:
    """Return canonical rate key strings for all reactions (for debugging)."""
    return [r.rate_key for r in self._reactions]

simulate(initial_conditions, t_span, t_eval=None, rtol=1e-08, atol=1e-12)

Integrate the ODE system forward in time and return the full trajectory.

Parameters

initial_conditions : dict[str, float] Initial concentrations; missing species default to 0. t_span : (t0, tf) Start and end times in seconds. t_eval : array-like, optional Times at which to store the solution. Defaults to 200 log-spaced points across t_span. rtol, atol : float Solver tolerances (passed to scipy Radau).

Returns

SimulationResult .times, .concentrations (dict species → 1-D array), .success. Use .final() for the last time-point dict or .at(t) for a specific time.

Source code in mantis/network.py

def simulate(
    self,
    initial_conditions: dict[str, float],
    t_span: tuple[float, float],
    t_eval=None,
    rtol: float = 1e-8,
    atol: float = 1e-12,
):
    """
    Integrate the ODE system forward in time and return the full trajectory.

    Parameters
    ----------
    initial_conditions : dict[str, float]
        Initial concentrations; missing species default to 0.
    t_span : (t0, tf)
        Start and end times in seconds.
    t_eval : array-like, optional
        Times at which to store the solution.  Defaults to 200 log-spaced
        points across t_span.
    rtol, atol : float
        Solver tolerances (passed to scipy Radau).

    Returns
    -------
    SimulationResult
        ``.times``, ``.concentrations`` (dict species → 1-D array),
        ``.success``.  Use ``.final()`` for the last time-point dict
        or ``.at(t)`` for a specific time.
    """
    from .analysis import simulate_ode
    return simulate_ode(
        self._reactions,
        self.species,
        self._rates,
        initial_conditions,
        t_span=t_span,
        t_eval=t_eval,
        chemostatted_values=self._chemostatted or None,
        rtol=rtol,
        atol=atol,
    )

stationary_distribution(initial_conditions, volume_L, *, initial_as='concentration', max_states=200000)

Exact stochastic stationary distribution (Anderson–Craciun–Kurtz 2010).

For a complex-balanced network the chemical master equation has a closed-form product-of-Poissons stationary distribution with means λ_i = N_A·V·c_i* (c* the complex-balanced equilibrium) — written down from structure alone, with no simulation. With conservation laws it is that product conditioned on the conserved totals (e.g. a Binomial for a single moiety).

Parameters

initial_conditions : dict[str, float] Fixes the compatibility class; interpreted per initial_as. volume_L : float Reaction volume in litres (sets the absolute molecule counts). initial_as : 'concentration' or 'count' max_states : int Cap on enumerated states for the conditioned (closed-network) case.

Returns

StationaryDistribution With .poisson_means(), .probability(state), .expected_counts(), .marginal(species) and .sample(...).

Raises

ValueError If the network is not complex-balanced (the closed form does not apply).

Source code in mantis/network.py

def stationary_distribution(
    self,
    initial_conditions: dict[str, float],
    volume_L: float,
    *,
    initial_as: str = "concentration",
    max_states: int = 200_000,
):
    """
    Exact stochastic stationary distribution (Anderson–Craciun–Kurtz 2010).

    For a complex-balanced network the chemical master equation has a closed-form
    product-of-Poissons stationary distribution with means ``λ_i = N_A·V·c_i*``
    (``c*`` the complex-balanced equilibrium) — written down *from structure
    alone, with no simulation*. With conservation laws it is that product
    conditioned on the conserved totals (e.g. a Binomial for a single moiety).

    Parameters
    ----------
    initial_conditions : dict[str, float]
        Fixes the compatibility class; interpreted per ``initial_as``.
    volume_L : float
        Reaction volume in litres (sets the absolute molecule counts).
    initial_as : 'concentration' or 'count'
    max_states : int
        Cap on enumerated states for the conditioned (closed-network) case.

    Returns
    -------
    StationaryDistribution
        With ``.poisson_means()``, ``.probability(state)``, ``.expected_counts()``,
        ``.marginal(species)`` and ``.sample(...)``.

    Raises
    ------
    ValueError
        If the network is not complex-balanced (the closed form does not apply).
    """
    from .stochastic_stationary import stationary_distribution
    return stationary_distribution(
        self._reactions, self.species, self._rates,
        self.deficiency, self.is_weakly_reversible,
        initial_conditions, volume_L,
        initial_as=initial_as,
        chemostatted_values=self._chemostatted or None,
        max_states=max_states,
    )

steady_states(initial_conditions, n_attempts=50, seed=None, t_end=10000.0)

Find steady states from the given initial conditions.

Uses ODE integration (Radau) as the primary strategy, which exactly preserves conservation laws, then falls back to multi-start least_squares for additional steady states. Duplicate solutions (relative L2 distance < 10%) and non-physical solutions (any concentration < −tol) are discarded. Results are sorted by residual ascending and filtered by relative residual against the best solution found to remove spurious stuck-at-IC artifacts.

Parameters

initial_conditions : dict[str, float] Initial concentrations; missing species default to 0. n_attempts : int Total number of solver attempts (including the integration from IC). seed : int, optional RNG seed for reproducibility. t_end : float ODE integration horizon in seconds (default 1e4). Increase for systems with slow reactions whose timescale τ = 1/(k·[X]) >> 1e4 s.

Returns

list[SteadyState] Each entry has .concentrations, .eigenvalues, .is_stable, .is_oscillatory, and .residual.

Source code in mantis/network.py

def steady_states(
    self,
    initial_conditions: dict[str, float],
    n_attempts: int = 50,
    seed: int | None = None,
    t_end: float = 1e4,
) -> list:
    """
    Find steady states from the given initial conditions.

    Uses ODE integration (Radau) as the primary strategy, which exactly
    preserves conservation laws, then falls back to multi-start least_squares
    for additional steady states. Duplicate solutions (relative L2 distance
    < 10%) and non-physical solutions (any concentration < −tol) are discarded.
    Results are sorted by residual ascending and filtered by relative residual
    against the best solution found to remove spurious stuck-at-IC artifacts.

    Parameters
    ----------
    initial_conditions : dict[str, float]
        Initial concentrations; missing species default to 0.
    n_attempts : int
        Total number of solver attempts (including the integration from IC).
    seed : int, optional
        RNG seed for reproducibility.
    t_end : float
        ODE integration horizon in seconds (default 1e4).  Increase for
        systems with slow reactions whose timescale τ = 1/(k·[X]) >> 1e4 s.

    Returns
    -------
    list[SteadyState]
        Each entry has ``.concentrations``, ``.eigenvalues``, ``.is_stable``,
        ``.is_oscillatory``, and ``.residual``.
    """
    from .analysis import find_steady_states
    return find_steady_states(
        self._reactions,
        self.species,
        self._rates,
        initial_conditions,
        n_attempts=n_attempts,
        seed=seed,
        chemostatted_values=self._chemostatted or None,
        t_end=t_end,
    )

stochastic_simulate(initial_conditions, t_span, volume_L, *, initial_as='concentration', max_events=1000000, seed=None)

Single-trajectory Gillespie SSA realization.

Use when molecule counts are low (∼ ≤ 10³) and the deterministic ODE gives the wrong answer — e.g., a CHA cascade at single-cell concentrations or stochastic switching in a bistable circuit.

Parameters

initial_conditions : species → initial count or concentration. t_span : (t0, tf) in seconds. volume_L : reaction volume in liters (e.g. 1e-4 = 100 µL). initial_as : 'concentration' (default) or 'count'. max_events : safety cap on reaction firings. seed : RNG seed for reproducibility.

Returns

StochasticResult with .times, .counts, .concentrations, .at(t), .final().

Source code in mantis/network.py

def stochastic_simulate(
    self,
    initial_conditions: dict[str, float],
    t_span: tuple[float, float],
    volume_L: float,
    *,
    initial_as: str = "concentration",
    max_events: int = 1_000_000,
    seed: int | None = None,
):
    """
    Single-trajectory Gillespie SSA realization.

    Use when molecule counts are low (∼ ≤ 10³) and the deterministic ODE
    gives the wrong answer — e.g., a CHA cascade at single-cell
    concentrations or stochastic switching in a bistable circuit.

    Parameters
    ----------
    initial_conditions : species → initial count or concentration.
    t_span             : (t0, tf) in seconds.
    volume_L           : reaction volume in liters (e.g. 1e-4 = 100 µL).
    initial_as         : 'concentration' (default) or 'count'.
    max_events         : safety cap on reaction firings.
    seed               : RNG seed for reproducibility.

    Returns
    -------
    StochasticResult with ``.times``, ``.counts``, ``.concentrations``,
    ``.at(t)``, ``.final()``.
    """
    from .analysis import gillespie_simulate
    return gillespie_simulate(
        self._reactions,
        self.species,
        self._rates,
        initial_conditions,
        t_span=t_span,
        volume_L=volume_L,
        initial_as=initial_as,
        max_events=max_events,
        seed=seed,
        chemostatted_values=self._chemostatted or None,
    )

tau_leap_simulate(initial_conditions, t_span, volume_L, *, initial_as='concentration', tau=None, epsilon=0.03, n_record=200, seed=None)

τ-leap stochastic simulation (approximate Gillespie).

Same interface as :meth:stochastic_simulate but fires reactions in Poisson-distributed bursts over each leap rather than one at a time. Roughly N× faster than direct SSA when populations are large (N is the mean number of firings per leap), at the cost of asymptotic accuracy as populations shrink.

Parameters

tau : optional fixed leap size (s); if None, adaptive τ per Cao 2006. epsilon : adaptive-τ tolerance (max fractional propensity change per leap). Smaller = more accurate, slower. n_record : number of evenly-spaced time points to record.

See :func:mantis.analysis.tau_leap_simulate for full details.

Source code in mantis/network.py

def tau_leap_simulate(
    self,
    initial_conditions: dict[str, float],
    t_span: tuple[float, float],
    volume_L: float,
    *,
    initial_as: str = "concentration",
    tau: float | None = None,
    epsilon: float = 0.03,
    n_record: int = 200,
    seed: int | None = None,
):
    """
    τ-leap stochastic simulation (approximate Gillespie).

    Same interface as :meth:`stochastic_simulate` but fires reactions in
    Poisson-distributed bursts over each leap rather than one at a time.
    Roughly N× faster than direct SSA when populations are large (N is
    the mean number of firings per leap), at the cost of asymptotic
    accuracy as populations shrink.

    Parameters
    ----------
    tau : optional fixed leap size (s); if None, adaptive τ per Cao 2006.
    epsilon : adaptive-τ tolerance (max fractional propensity change per
              leap).  Smaller = more accurate, slower.
    n_record : number of evenly-spaced time points to record.

    See :func:`mantis.analysis.tau_leap_simulate` for full details.
    """
    from .analysis import tau_leap_simulate
    return tau_leap_simulate(
        self._reactions,
        self.species,
        self._rates,
        initial_conditions,
        t_span=t_span,
        volume_L=volume_L,
        initial_as=initial_as,
        tau=tau,
        epsilon=epsilon,
        n_record=n_record,
        seed=seed,
        chemostatted_values=self._chemostatted or None,
    )

to_sbml(path=None)

Export this network to SBML (Level 3 Version 2).

Requires the optional python-libsbml package (pip install mantis-delta[sbml]).

Parameters

path : str, optional If given, the SBML document is written to this file path. The XML string is returned in all cases.

Returns

str The SBML document as an XML string.

Source code in mantis/network.py

def to_sbml(self, path: str | None = None) -> str:
    """
    Export this network to SBML (Level 3 Version 2).

    Requires the optional ``python-libsbml`` package
    (``pip install mantis-delta[sbml]``).

    Parameters
    ----------
    path : str, optional
        If given, the SBML document is written to this file path.  The XML
        string is returned in all cases.

    Returns
    -------
    str
        The SBML document as an XML string.
    """
    from .sbml import network_to_sbml, write_sbml
    if path is not None:
        write_sbml(self, path)
    return network_to_sbml(self)

mantis.analysis

Numerical ODE simulation, steady-state finding, stability analysis, and bifurcation scanning.

SimulationResult dataclass

Source code in mantis/analysis.py

@dataclass
class SimulationResult:
    times: np.ndarray
    concentrations: dict[str, np.ndarray]
    success: bool

    def at(self, t: float) -> dict[str, float]:
        """Return interpolated concentrations at time t."""
        idx = int(np.searchsorted(self.times, t, side="right")) - 1
        idx = max(0, min(idx, len(self.times) - 1))
        return {sp: float(arr[idx]) for sp, arr in self.concentrations.items()}

    def final(self) -> dict[str, float]:
        """Return concentrations at the last time point."""
        return {sp: float(arr[-1]) for sp, arr in self.concentrations.items()}

at(t)

Return interpolated concentrations at time t.

Source code in mantis/analysis.py

def at(self, t: float) -> dict[str, float]:
    """Return interpolated concentrations at time t."""
    idx = int(np.searchsorted(self.times, t, side="right")) - 1
    idx = max(0, min(idx, len(self.times) - 1))
    return {sp: float(arr[idx]) for sp, arr in self.concentrations.items()}

final()

Return concentrations at the last time point.

Source code in mantis/analysis.py

def final(self) -> dict[str, float]:
    """Return concentrations at the last time point."""
    return {sp: float(arr[-1]) for sp, arr in self.concentrations.items()}

StochasticResult dataclass

Trajectory from a single Gillespie SSA realization.

Attributes

times : 1-D ndarray of reaction times (length n_events + 1) counts : dict species → 1-D ndarray of integer molecule counts concentrations : dict species → 1-D ndarray of concentrations (M) (counts / (volume * Avogadro)) n_events : number of reaction firings recorded success : True if integration finished before exhausting events volume_L : reaction volume in liters

Source code in mantis/analysis.py

@dataclass
class StochasticResult:
    """Trajectory from a single Gillespie SSA realization.

    Attributes
    ----------
    times          : 1-D ndarray of reaction times (length n_events + 1)
    counts         : dict species → 1-D ndarray of integer molecule counts
    concentrations : dict species → 1-D ndarray of concentrations (M)
                     (counts / (volume * Avogadro))
    n_events       : number of reaction firings recorded
    success        : True if integration finished before exhausting events
    volume_L       : reaction volume in liters
    """
    times: np.ndarray
    counts: dict[str, np.ndarray]
    concentrations: dict[str, np.ndarray]
    n_events: int
    success: bool
    volume_L: float

    def at(self, t: float) -> dict[str, float]:
        """Step-function lookup of concentrations at time t."""
        idx = int(np.searchsorted(self.times, t, side="right")) - 1
        idx = max(0, min(idx, len(self.times) - 1))
        return {sp: float(arr[idx]) for sp, arr in self.concentrations.items()}

    def final(self) -> dict[str, float]:
        return {sp: float(arr[-1]) for sp, arr in self.concentrations.items()}

at(t)

Step-function lookup of concentrations at time t.

Source code in mantis/analysis.py

def at(self, t: float) -> dict[str, float]:
    """Step-function lookup of concentrations at time t."""
    idx = int(np.searchsorted(self.times, t, side="right")) - 1
    idx = max(0, min(idx, len(self.times) - 1))
    return {sp: float(arr[idx]) for sp, arr in self.concentrations.items()}

build_ode_function(reactions, species, rate_values, chemostatted=None)

Return a callable f(t, y) → dydt using pure numpy. Pre-computes N matrix and ordered rate array for speed.

Chemostatted species are excluded from reactants_info (their contribution has already been folded into rate_values via fold_chemostatted_into_rates).

Source code in mantis/analysis.py

def build_ode_function(
    reactions: list[Reaction],
    species: list[str],
    rate_values: dict[str, float],
    chemostatted: set[str] | None = None,
) -> Any:
    """
    Return a callable f(t, y) → dydt using pure numpy.
    Pre-computes N matrix and ordered rate array for speed.

    Chemostatted species are excluded from reactants_info (their contribution
    has already been folded into rate_values via fold_chemostatted_into_rates).
    """
    chem = chemostatted or set()
    sp_idx = {s: i for i, s in enumerate(species)}
    N = build_stoichiometry_matrix(reactions, species)
    rates_arr = np.array([rate_values.get(r.rate_key, 0.0) for r in reactions])

    # Pre-compute reactant info: list of (species_index, coeff) per reaction
    # Skip chemostatted species — their concentrations are already in rates_arr
    reactants_info = [
        [(sp_idx[name], coeff) for name, coeff in rxn.reactants if name not in chem]
        for rxn in reactions
    ]

    def f(t, y):
        y = np.maximum(y, 0.0)
        fluxes = np.empty(len(reactions))
        for j, (rate, ri) in enumerate(zip(rates_arr, reactants_info)):
            flux = rate
            for idx, coeff in ri:
                flux *= y[idx] ** coeff
            fluxes[j] = flux
        return N @ fluxes

    return f

classify_steady_state(eigenvalues, tol_zero=0.0001)

Returns (is_stable, is_oscillatory).

is_stable — True iff all significant eigenvalues have Re < 0. is_oscillatory — True iff any significant eigenvalue has a non-trivial imaginary part (|Im| > 1e-10 * |λ|). This covers both stable spirals (Re<0) and unstable foci/spirals (Re>0), i.e. any fixed point near a Hopf bifurcation.

Filters near-zero eigenvalues (from conservation law dimensions) before classifying.

Source code in mantis/analysis.py

def classify_steady_state(
    eigenvalues: np.ndarray,
    tol_zero: float = 1e-4,
) -> tuple[bool, bool]:
    """
    Returns (is_stable, is_oscillatory).

    ``is_stable``     — True iff all significant eigenvalues have Re < 0.
    ``is_oscillatory`` — True iff any significant eigenvalue has a non-trivial
                         imaginary part (|Im| > 1e-10 * |λ|).  This covers both
                         stable spirals (Re<0) *and* unstable foci/spirals (Re>0),
                         i.e. any fixed point near a Hopf bifurcation.

    Filters near-zero eigenvalues (from conservation law dimensions) before classifying.
    """
    if len(eigenvalues) == 0:
        return True, False
    max_abs = np.max(np.abs(eigenvalues))
    if max_abs < 1e-30:
        return True, False
    # Filter eigenvalues that are effectively zero relative to the largest
    significant = eigenvalues[np.abs(eigenvalues) > tol_zero * max_abs]
    if len(significant) == 0:
        return True, False
    is_stable = bool(np.all(np.real(significant) < 0))
    # Oscillatory = complex eigenvalues regardless of stability
    # (captures both stable spirals and unstable foci / Hopf-born limit cycles)
    is_oscillatory = bool(
        np.any(np.abs(np.imag(significant)) > 1e-10 * np.abs(significant))
    )
    return is_stable, is_oscillatory

find_steady_states(reactions, species, rate_values, initial_conditions, n_attempts=50, tol=1e-08, seed=None, chemostatted_values=None, t_end=10000.0)

Steady-state finder using three strategies: 1. Direct least_squares from user IC (finds both stable and unstable fixed points). 2. ODE integration from user IC to t_end (conserves CLs; reaches stable attractors). 3. Multi-start with scale-aware random ICs → ODE then least_squares fallback.

Chemostatted species are folded into effective rate constants so the ODE system only tracks dynamic species.

Parameters

t_end : float Integration horizon for ODE-based strategies (seconds). The default (1e4 s) is intentionally short so that oscillatory / limit-cycle systems do not hang; the algebraic least_squares fallback handles those cases. For systems with very slow reactions (e.g. leakage pathways with τ >> 1e4 s) increase this value so that ODE integration reaches the true attractor. A safe upper bound is 10–100× the slowest relevant timescale (τ = 1 / (k_slow × [reactant])).

Returns a list of SteadyState objects, sorted by true ODE residual (lowest first), filtered to remove duplicate states and high-residual algebraic artifacts.

Source code in mantis/analysis.py

def find_steady_states(
    reactions: list[Reaction],
    species: list[str],
    rate_values: dict[str, float],
    initial_conditions: dict[str, float],
    n_attempts: int = 50,
    tol: float = 1e-8,
    seed: int | None = None,
    chemostatted_values: dict[str, float] | None = None,
    t_end: float = 1e4,
) -> list[SteadyState]:
    """
    Steady-state finder using three strategies:
    1. Direct least_squares from user IC (finds both stable and unstable fixed points).
    2. ODE integration from user IC to t_end (conserves CLs; reaches stable attractors).
    3. Multi-start with scale-aware random ICs → ODE then least_squares fallback.

    Chemostatted species are folded into effective rate constants so the ODE
    system only tracks dynamic species.

    Parameters
    ----------
    t_end : float
        Integration horizon for ODE-based strategies (seconds).  The default
        (1e4 s) is intentionally short so that oscillatory / limit-cycle
        systems do not hang; the algebraic least_squares fallback handles those
        cases.  For systems with very slow reactions (e.g. leakage pathways
        with τ >> 1e4 s) increase this value so that ODE integration reaches
        the true attractor.  A safe upper bound is 10–100× the slowest
        relevant timescale (τ = 1 / (k_slow × [reactant])).

    Returns a list of SteadyState objects, sorted by true ODE residual (lowest first),
    filtered to remove duplicate states and high-residual algebraic artifacts.
    """
    from .stoichiometry import fold_chemostatted_into_rates
    chem_vals = chemostatted_values or {}
    chem_keys = set(chem_vals.keys())
    eff_rates = fold_chemostatted_into_rates(reactions, rate_values, chem_vals) if chem_vals else rate_values

    n_sp = len(species)
    rng = np.random.default_rng(seed)

    cl_vecs = _conservation_law_vectors(reactions, species)
    y_ref = np.array([initial_conditions.get(s, 0.0) for s in species])
    cl_totals = [float(v @ y_ref) for v in cl_vecs]

    f_ode = build_ode_function(reactions, species, eff_rates, chem_keys)
    full_jac = _full_jacobian_fn(reactions, species, eff_rates, chem_keys)
    system_fn = _build_augmented_system(reactions, species, eff_rates, cl_vecs, cl_totals, chem_keys)

    collected: list[SteadyState] = []

    def _try_add(y_sol: np.ndarray, _residual_hint: float) -> None:
        y_sol = np.maximum(y_sol, 0.0)
        # Recompute the true ODE residual |dy/dt|.  We store it for post-filtering
        # (relative to the best-found solution) rather than applying a fixed absolute
        # threshold, which fails on stiff systems where all flux magnitudes are tiny.
        abs_res = float(np.linalg.norm(f_ode(0, y_sol)))
        if np.any(y_sol < -tol):
            return
        for existing in collected:
            y_ex = np.array([existing.concentrations[s] for s in species])
            if np.linalg.norm(y_sol - y_ex) / (np.linalg.norm(y_ex) + 1e-30) < 0.10:
                return
        J_num = full_jac(y_sol)
        eigs = eigvals(J_num)
        is_stable, is_osc = classify_steady_state(eigs)
        ss = SteadyState(
            concentrations=dict(zip(species, y_sol.tolist())),
            eigenvalues=eigs,
            is_stable=is_stable,
            is_oscillatory=is_osc,
            residual=abs_res,
        )
        collected.append(ss)

    def _try_least_squares(y0: np.ndarray) -> None:
        """Try algebraic least_squares from y0; can find unstable fixed points."""
        try:
            result = least_squares(
                system_fn,
                y0,
                bounds=(0.0, np.inf),
                method="trf",
                ftol=tol,
                xtol=tol,
                gtol=tol,
                max_nfev=5000,
            )
            _try_add(result.x, float(np.linalg.norm(result.fun)))
        except Exception:
            pass

    has_cl = len(cl_vecs) > 0

    if has_cl:
        # ── Closed / semi-closed system (has conservation laws) ──────────────
        # Strategy 1a: ODE integration from user IC — respects CLs, reaches the
        # physically relevant attractor (unique within each conservation class).
        y_final, res = _integrate_to_ss(f_ode, y_ref.copy(), t_end=t_end)
        if y_final is not None:
            _try_add(y_final, res)

        # Strategy 1b: algebraic solve from user IC — finds any fixed point
        # including those that are unstable (integration diverges from them).
        _try_least_squares(y_ref.copy())

        # Strategy 2: multi-start on the constraint manifold
        scale = max(float(np.max(np.abs(y_ref))), 1.0)
        for _ in range(n_attempts - 2):
            y0 = _make_feasible_initial(cl_vecs, cl_totals, n_sp, rng)
            if y0 is None:
                continue
            y_ivp, res_ivp = _integrate_to_ss(f_ode, y0, t_end=t_end)
            if y_ivp is not None and res_ivp < 1e-4:
                _try_add(y_ivp, res_ivp)
                continue
            _try_least_squares(y0)

    else:
        # ── Open / chemostatted system (no conservation laws) ─────────────────
        # ODE integration is useless here: limit-cycle attractors make it hang,
        # and there is no conservation manifold to project onto.  Use algebraic
        # least_squares exclusively.

        # Strategy 1: algebraic solve from user IC (finds both stable & unstable FPs)
        _try_least_squares(y_ref.copy())

        # Strategy 2: scale-aware multi-start algebraic solve
        scale = max(float(np.max(np.abs(y_ref))), 1.0)
        for _ in range(n_attempts - 1):
            y0 = np.abs(rng.normal(loc=scale, scale=scale * 0.5, size=n_sp))
            y0 = np.maximum(y0, 1e-12)
            _try_least_squares(y0)

    collected.sort(key=lambda s: s.residual)
    # Post-filter: drop candidates whose ODE residual is more than 1e4× worse than
    # the best (lowest-residual) solution.  This correctly handles stiff systems
    # where all fluxes are tiny in absolute terms — a relative comparison is the
    # only meaningful discriminant between "converged" and "stuck-at-IC" solutions.
    if collected:
        best_res = collected[0].residual
        collected = [s for s in collected if s.residual <= best_res * 1e4 + 1e-30]
    return collected

gillespie_simulate(reactions, species, rate_values, initial_conditions, t_span, volume_L, *, initial_as='concentration', max_events=1000000, seed=None, chemostatted_values=None)

Single-trajectory Gillespie direct-method SSA.

Parameters

reactions, species, rate_values Same shape as :func:simulate_ode — deterministic mass-action rates. initial_conditions Either molecule counts (set initial_as='count') or molar concentrations (default). t_span : (t0, tf) Simulation interval. volume_L : float Reaction volume in liters. Required because stochastic propensities depend on absolute molecule counts. For 100 µL physiological volume, pass 1e-4. initial_as : 'concentration' or 'count' max_events : safety cap on reaction firings. seed : optional RNG seed for reproducibility. chemostatted_values : species kept at constant concentration; folded into the propensities of reactions that consume them.

Returns

StochasticResult Full trajectory of times, counts and concentrations.

Source code in mantis/analysis.py

def gillespie_simulate(
    reactions: list[Reaction],
    species: list[str],
    rate_values: dict[str, float],
    initial_conditions: dict[str, float],
    t_span: tuple[float, float],
    volume_L: float,
    *,
    initial_as: str = "concentration",
    max_events: int = 1_000_000,
    seed: int | None = None,
    chemostatted_values: dict[str, float] | None = None,
) -> StochasticResult:
    """
    Single-trajectory Gillespie direct-method SSA.

    Parameters
    ----------
    reactions, species, rate_values
        Same shape as :func:`simulate_ode` — deterministic mass-action rates.
    initial_conditions
        Either molecule counts (set ``initial_as='count'``) or molar
        concentrations (default).
    t_span : (t0, tf)
        Simulation interval.
    volume_L : float
        Reaction volume in liters.  Required because stochastic propensities
        depend on absolute molecule counts.  For 100 µL physiological volume,
        pass ``1e-4``.
    initial_as : 'concentration' or 'count'
    max_events : safety cap on reaction firings.
    seed : optional RNG seed for reproducibility.
    chemostatted_values : species kept at constant concentration; folded into
        the propensities of reactions that consume them.

    Returns
    -------
    StochasticResult
        Full trajectory of times, counts and concentrations.
    """
    from .stoichiometry import fold_chemostatted_into_rates
    from ._kernels import build_kernel_arrays, normalize_seed, ssa_kernel

    chem_vals = chemostatted_values or {}
    eff_rates = (
        fold_chemostatted_into_rates(reactions, rate_values, chem_vals)
        if chem_vals else rate_values
    )
    chem = set(chem_vals.keys())

    sp_idx = {s: i for i, s in enumerate(species)}
    n_sp = len(species)
    NA_V = AVOGADRO * volume_L

    # Initial counts
    n0 = np.zeros(n_sp, dtype=np.int64)
    for s, val in initial_conditions.items():
        if s not in sp_idx:
            continue
        if initial_as == "count":
            n0[sp_idx[s]] = int(round(val))
        else:
            n0[sp_idx[s]] = int(round(val * NA_V))

    c_eff, react_sp, react_co, change = build_kernel_arrays(
        reactions, species, eff_rates, chem, NA_V
    )

    # The direct-method SSA loop runs in a numba-JIT kernel when available,
    # else as an equivalent pure-Python loop (see mantis._kernels).
    times_arr, history_arr, n_events, success = ssa_kernel(
        n0, c_eff, react_sp, react_co, change,
        float(t_span[0]), float(t_span[1]), int(max_events), normalize_seed(seed),
    )

    counts_dict = {sp: history_arr[:, i] for i, sp in enumerate(species)}
    conc_dict = {sp: history_arr[:, i] / NA_V for i, sp in enumerate(species)}

    return StochasticResult(
        times=times_arr,
        counts=counts_dict,
        concentrations=conc_dict,
        n_events=int(n_events),
        success=bool(success),
        volume_L=volume_L,
    )

scan_bifurcation(reactions, species, base_rate_values, parameter, param_range, n_points, initial_conditions, n_attempts=20, chemostatted_values=None, t_end=10000.0)

Vary one rate constant over a log-spaced range and collect steady states.

Parameters

t_end : float ODE integration horizon passed to find_steady_states at each point. See that function's t_end documentation for guidance on choosing a value appropriate for the slowest reaction timescale in your network.

Source code in mantis/analysis.py

def scan_bifurcation(
    reactions: list[Reaction],
    species: list[str],
    base_rate_values: dict[str, float],
    parameter: str,
    param_range: tuple[float, float],
    n_points: int,
    initial_conditions: dict[str, float],
    n_attempts: int = 20,
    chemostatted_values: dict[str, float] | None = None,
    t_end: float = 1e4,
) -> BifurcationResult:
    """
    Vary one rate constant over a log-spaced range and collect steady states.

    Parameters
    ----------
    t_end : float
        ODE integration horizon passed to ``find_steady_states`` at each point.
        See that function's ``t_end`` documentation for guidance on choosing a
        value appropriate for the slowest reaction timescale in your network.
    """
    from .parsing import normalize_rate_key
    try:
        norm_param = normalize_rate_key(parameter)
    except ValueError:
        norm_param = parameter

    param_values = np.logspace(
        np.log10(param_range[0]), np.log10(param_range[1]), n_points
    )
    all_ss: list[list[SteadyState]] = []
    for pval in param_values:
        rates = dict(base_rate_values)
        rates[norm_param] = pval
        ss_list = find_steady_states(
            reactions, species, rates, initial_conditions,
            n_attempts=n_attempts, seed=42,
            chemostatted_values=chemostatted_values,
            t_end=t_end,
        )
        all_ss.append(ss_list)

    return BifurcationResult(
        parameter_name=norm_param,
        parameter_values=param_values,
        steady_states=all_ss,
    )

simulate_ode(reactions, species, rate_values, initial_conditions, t_span, t_eval=None, chemostatted_values=None, rtol=1e-08, atol=1e-12)

Integrate the ODE system forward in time and return the full trajectory.

Parameters

t_span : (t0, tf) Start and end times in seconds. t_eval : array-like, optional Times at which to store the solution. Defaults to 200 log-spaced points across t_span (or linear if t_span[0] == 0 and t_span[1] <= 0). chemostatted_values : dict, optional Fixed concentrations folded into rate constants (same semantics as find_steady_states).

Returns

SimulationResult .times (1-D array), .concentrations (dict species → 1-D array), .success (bool). Use .final() to get the last time-point dict or .at(t) for a specific time.

Source code in mantis/analysis.py

def simulate_ode(
    reactions: list[Reaction],
    species: list[str],
    rate_values: dict[str, float],
    initial_conditions: dict[str, float],
    t_span: tuple[float, float],
    t_eval: np.ndarray | None = None,
    chemostatted_values: dict[str, float] | None = None,
    rtol: float = 1e-8,
    atol: float = 1e-12,
) -> SimulationResult:
    """
    Integrate the ODE system forward in time and return the full trajectory.

    Parameters
    ----------
    t_span : (t0, tf)
        Start and end times in seconds.
    t_eval : array-like, optional
        Times at which to store the solution.  Defaults to 200 log-spaced
        points across t_span (or linear if t_span[0] == 0 and t_span[1] <= 0).
    chemostatted_values : dict, optional
        Fixed concentrations folded into rate constants (same semantics as
        ``find_steady_states``).

    Returns
    -------
    SimulationResult
        ``.times`` (1-D array), ``.concentrations`` (dict species → 1-D array),
        ``.success`` (bool).  Use ``.final()`` to get the last time-point dict
        or ``.at(t)`` for a specific time.
    """
    from scipy.integrate import solve_ivp
    from .stoichiometry import fold_chemostatted_into_rates

    chem_vals = chemostatted_values or {}
    eff_rates = fold_chemostatted_into_rates(reactions, rate_values, chem_vals) if chem_vals else rate_values
    chem_keys = set(chem_vals.keys())

    f_ode = build_ode_function(reactions, species, eff_rates, chem_keys)
    y0 = np.array([initial_conditions.get(s, 0.0) for s in species])

    t0, tf = float(t_span[0]), float(t_span[1])
    if t_eval is None:
        n_pts = 200
        if t0 <= 0 and tf > 0:
            t_eval = np.logspace(np.log10(max(tf * 1e-6, 1e-6)), np.log10(tf), n_pts - 1)
            t_eval[-1] = tf
            t_eval = np.concatenate([[t0], t_eval])
        else:
            t_eval = np.linspace(t0, tf, n_pts)

    try:
        sol = solve_ivp(
            f_ode,
            [t0, tf],
            y0,
            method="Radau",
            t_eval=t_eval,
            rtol=rtol,
            atol=atol,
            dense_output=False,
        )
        times = sol.t
        conc = {sp: np.maximum(sol.y[i], 0.0) for i, sp in enumerate(species)}
        if not sol.success:
            print("solve_ivp failed:", getattr(sol, 'message', 'No message'))
        return SimulationResult(times=times, concentrations=conc, success=sol.success)
    except Exception as e:
        import traceback
        traceback.print_exc()
        times = np.array([t0, tf])
        conc = {sp: np.full(2, y0[i]) for i, sp in enumerate(species)}
        return SimulationResult(times=times, concentrations=conc, success=False)

tau_leap_simulate(reactions, species, rate_values, initial_conditions, t_span, volume_L, *, initial_as='concentration', tau=None, epsilon=0.03, n_record=200, seed=None, chemostatted_values=None)

Approximate Gillespie via τ-leap (Gillespie 2001; Cao, Gillespie, Petzold 2006).

Instead of one reaction per step, fire all reactions in parallel over a timestep τ, sampling each firing count from Poisson(a_j · τ). Much faster than direct SSA when populations are large and propensities change slowly, at the cost of accuracy: τ-leap is exact only in the limit of frequent reactions per leap (large molecule counts).

Parameters

reactions, species, rate_values, initial_conditions, t_span, volume_L, initial_as, seed, chemostatted_values Same as :func:gillespie_simulate. tau : optional fixed leap size (seconds). If None, an adaptive τ is chosen each step bounded by epsilon (Cao 2006 algorithm). epsilon : adaptive-τ tolerance. Controls the maximum allowed fractional change in any propensity per leap; smaller = more accurate but slower. Ignored when tau is given. n_record : approximate number of evenly-spaced time points to record in the returned trajectory (the simulator's internal step is independent of this — recording is via linear interpolation of step boundaries).

Returns

StochasticResult Same shape as :func:gillespie_simulate.

Notes

To preserve non-negativity in stiff systems, this implementation falls back to a single direct-SSA step when the chosen τ would drive any species count below zero (a simple "leap rejection" — not the binomial Tian/Burrage refinement, but adequate for most kinetic-design contexts).

Source code in mantis/analysis.py

def tau_leap_simulate(
    reactions: list[Reaction],
    species: list[str],
    rate_values: dict[str, float],
    initial_conditions: dict[str, float],
    t_span: tuple[float, float],
    volume_L: float,
    *,
    initial_as: str = "concentration",
    tau: float | None = None,
    epsilon: float = 0.03,
    n_record: int = 200,
    seed: int | None = None,
    chemostatted_values: dict[str, float] | None = None,
) -> StochasticResult:
    """
    Approximate Gillespie via τ-leap (Gillespie 2001; Cao, Gillespie, Petzold 2006).

    Instead of one reaction per step, fire all reactions in parallel over a
    timestep τ, sampling each firing count from Poisson(a_j · τ).  Much faster
    than direct SSA when populations are large and propensities change slowly,
    at the cost of accuracy: τ-leap is exact only in the limit of frequent
    reactions per leap (large molecule counts).

    Parameters
    ----------
    reactions, species, rate_values, initial_conditions, t_span, volume_L,
    initial_as, seed, chemostatted_values
        Same as :func:`gillespie_simulate`.
    tau : optional fixed leap size (seconds).  If None, an adaptive τ is
        chosen each step bounded by ``epsilon`` (Cao 2006 algorithm).
    epsilon : adaptive-τ tolerance.  Controls the maximum allowed fractional
        change in any propensity per leap; smaller = more accurate but slower.
        Ignored when ``tau`` is given.
    n_record : approximate number of evenly-spaced time points to record in
        the returned trajectory (the simulator's internal step is independent
        of this — recording is via linear interpolation of step boundaries).

    Returns
    -------
    StochasticResult
        Same shape as :func:`gillespie_simulate`.

    Notes
    -----
    To preserve non-negativity in stiff systems, this implementation falls
    back to a single direct-SSA step when the chosen τ would drive any
    species count below zero (a simple "leap rejection" — not the binomial
    Tian/Burrage refinement, but adequate for most kinetic-design contexts).
    """
    from .stoichiometry import fold_chemostatted_into_rates
    from ._kernels import build_kernel_arrays, normalize_seed, tau_leap_kernel

    chem_vals = chemostatted_values or {}
    eff_rates = (
        fold_chemostatted_into_rates(reactions, rate_values, chem_vals)
        if chem_vals else rate_values
    )
    chem = set(chem_vals.keys())

    sp_idx = {s: i for i, s in enumerate(species)}
    n_sp = len(species)
    NA_V = AVOGADRO * volume_L

    n0 = np.zeros(n_sp, dtype=np.int64)
    for s, val in initial_conditions.items():
        if s not in sp_idx:
            continue
        if initial_as == "count":
            n0[sp_idx[s]] = int(round(val))
        else:
            n0[sp_idx[s]] = int(round(val * NA_V))

    c_eff, react_sp, react_co, change_matrix = build_kernel_arrays(
        reactions, species, eff_rates, chem, NA_V
    )

    t0 = float(t_span[0])
    tf = float(t_span[1])
    record_times = np.linspace(t0, tf, max(2, n_record))

    # Adaptive (fixed_tau<=0) or fixed-step tau-leap, JIT-compiled when numba is
    # present, otherwise an equivalent pure-Python loop (see mantis._kernels).
    record_states, success = tau_leap_kernel(
        n0, c_eff, react_sp, react_co, change_matrix, t0, tf,
        float(tau) if tau is not None else -1.0, float(epsilon),
        record_times, normalize_seed(seed),
    )

    counts_dict = {sp: record_states[:, i] for i, sp in enumerate(species)}
    conc_dict = {sp: record_states[:, i] / NA_V for i, sp in enumerate(species)}

    return StochasticResult(
        times=record_times,
        counts=counts_dict,
        concentrations=conc_dict,
        n_events=int(len(record_times)),
        success=success,
        volume_L=volume_L,
    )

Exhaustive steady states

mantis.steady_states

Exhaustive steady-state enumeration for mass-action networks.

The hybrid solver in :mod:mantis.analysis (ODE integration + multi-start least_squares) can silently miss steady states — in particular unstable fixed points that forward integration flows away from. At mass-action equilibrium the positive steady states are the positive real roots of a polynomial system, so they can instead be found exhaustively with completeness guarantees:

• Default pure-Python backend (backend="sympy"). The reaction-rate polynomials are reduced against the conservation laws, a lexicographic Gröbner basis triangularises the system, every complex root is recovered by back-substitution (numpy.roots) and refined by a Newton polish, and the real non-negative roots are returned. For a zero-dimensional ideal this finds all solutions (their count is bounded by the Bézout number), so it certifies the full steady-state set rather than sampling it.
• Optional homotopy backend (backend="phcpy"). Delegates to PHCpack via :mod:phcpy (Verschelde) — numerical polynomial homotopy continuation that scales to larger / stiffer systems. Requires pip install mantis-delta[homotopy].

The classic demonstration is the Schlögl model: this enumerator returns all three positive steady states (including the unstable middle branch that integration never reaches), whereas the multi-start solver typically returns only the two stable ones.

all_steady_states(reactions, species, rate_values, initial_conditions, chemostatted_values=None, backend='auto', scale=None, residual_tol=1e-08, imag_tol=1e-06)

Exhaustively enumerate the positive steady states of a mass-action network.

Returns a list of :class:~mantis.analysis.SteadyState objects (with eigenvalues and stability classification), sorted by residual. See the module docstring for the algorithm and the meaning of backend.

Source code in mantis/steady_states.py

def all_steady_states(
    reactions: list[Reaction],
    species: list[str],
    rate_values: dict[str, float],
    initial_conditions: dict[str, float],
    chemostatted_values: dict[str, float] | None = None,
    backend: str = "auto",
    scale: float | None = None,
    residual_tol: float = 1e-8,
    imag_tol: float = 1e-6,
):
    """Exhaustively enumerate the positive steady states of a mass-action network.

    Returns a list of :class:`~mantis.analysis.SteadyState` objects (with
    eigenvalues and stability classification), sorted by residual.  See the module
    docstring for the algorithm and the meaning of ``backend``.
    """
    from scipy.linalg import eigvals
    from .analysis import SteadyState, _full_jacobian_fn, classify_steady_state, build_ode_function
    from .stoichiometry import fold_chemostatted_into_rates

    chem_vals = chemostatted_values or {}
    chem_keys = set(chem_vals.keys())
    eff_rates = (
        fold_chemostatted_into_rates(reactions, rate_values, chem_vals)
        if chem_vals else rate_values
    )

    if backend == "phcpy":
        roots = _phcpy_all_roots(reactions, species, rate_values, initial_conditions,
                                 chem_vals, scale, residual_tol, imag_tol)
    elif backend in ("auto", "sympy"):
        try:
            roots = _sympy_all_roots(reactions, species, rate_values, initial_conditions,
                                     chem_vals, scale, residual_tol, imag_tol)
        except Exception:
            if backend == "auto":
                try:
                    roots = _phcpy_all_roots(reactions, species, rate_values,
                                             initial_conditions, chem_vals, scale,
                                             residual_tol, imag_tol)
                except Exception:
                    raise
            else:
                raise
    else:
        raise ValueError(f"Unknown backend {backend!r}; use 'auto', 'sympy', or 'phcpy'.")

    f_ode = build_ode_function(reactions, species, eff_rates, chem_keys)
    full_jac = _full_jacobian_fn(reactions, species, eff_rates, chem_keys)

    states: list[SteadyState] = []
    for y in roots:
        eigs = eigvals(full_jac(y))
        is_stable, is_osc = classify_steady_state(eigs)
        states.append(SteadyState(
            concentrations=dict(zip(species, y.tolist())),
            eigenvalues=eigs,
            is_stable=is_stable,
            is_oscillatory=is_osc,
            residual=float(np.linalg.norm(f_ode(0, y))),
        ))
    states.sort(key=lambda s: s.residual)
    return states

Global stability (Horn–Jackson Lyapunov)

mantis.stability

Global asymptotic stability certification via the Horn–Jackson Lyapunov function.

The structural theorems already in mantis certify local stability through the sign of the Jacobian eigenvalues. For complex-balanced networks a much stronger, global statement is available (Horn & Jackson 1972; Feinberg's Deficiency Zero Theorem): the pseudo-Helmholtz function

V(c) = Σ_i [ c_i (ln(c_i / c_i*) − 1) + c_i* ]

is a strict Lyapunov function on every positive stoichiometric compatibility class, where c* is the (unique) positive complex-balanced equilibrium of that class. Indeed ∂V/∂c_i = ln(c_i/c_i*), so

dV/dt = Σ_i ln(c_i/c_i*) · ċ_i ≤ 0,

with equality only at c = c*. Hence c* is globally asymptotically stable within the relative interior of its compatibility class — strictly stronger than the eigenvalue check.

Every weakly reversible deficiency-zero network is complex-balanced for all rate constants (Horn–Jackson), so the certificate is purely structural there; for higher-deficiency weakly-reversible networks complex balance is rate-dependent and is verified numerically at the equilibrium.

StabilityCertificate dataclass

Outcome of a global-stability analysis.

Attributes

is_complex_balanced : bool Whether the network is complex-balanced (structurally for δ=0 weakly reversible, else numerically at the equilibrium). globally_stable : bool True when a Horn–Jackson Lyapunov certificate establishes global asymptotic stability within the compatibility class. equilibrium : dict[str, float] | None The positive complex-balanced equilibrium c* (per the supplied class). lyapunov_function : sympy.Expr | None The symbolic pseudo-Helmholtz function V(c). lyapunov_derivative : sympy.Expr | None The symbolic dV/dt = Σ ln(c_i/c_i*)·f_i(c). n_samples_checked : int Number of random points in the compatibility class at which dV/dt ≤ 0 was numerically verified. max_dVdt : float Largest dV/dt observed over the samples (should be ≤ 0). rationale : str Human-readable explanation of the verdict.

Source code in mantis/stability.py

@dataclass
class StabilityCertificate:
    """Outcome of a global-stability analysis.

    Attributes
    ----------
    is_complex_balanced : bool
        Whether the network is complex-balanced (structurally for δ=0 weakly
        reversible, else numerically at the equilibrium).
    globally_stable : bool
        True when a Horn–Jackson Lyapunov certificate establishes global
        asymptotic stability within the compatibility class.
    equilibrium : dict[str, float] | None
        The positive complex-balanced equilibrium ``c*`` (per the supplied class).
    lyapunov_function : sympy.Expr | None
        The symbolic pseudo-Helmholtz function ``V(c)``.
    lyapunov_derivative : sympy.Expr | None
        The symbolic ``dV/dt = Σ ln(c_i/c_i*)·f_i(c)``.
    n_samples_checked : int
        Number of random points in the compatibility class at which ``dV/dt ≤ 0``
        was numerically verified.
    max_dVdt : float
        Largest ``dV/dt`` observed over the samples (should be ≤ 0).
    rationale : str
        Human-readable explanation of the verdict.
    """
    is_complex_balanced: bool
    globally_stable: bool
    equilibrium: dict[str, float] | None = None
    lyapunov_function: "sympy.Expr | None" = None
    lyapunov_derivative: "sympy.Expr | None" = None
    n_samples_checked: int = 0
    max_dVdt: float = 0.0
    rationale: str = ""
    summary_lines: list[str] = field(default_factory=list)

    def __str__(self) -> str:
        return "\n".join(self.summary_lines) if self.summary_lines else self.rationale

certify_global_stability(reactions, species, rate_values, deficiency, is_weakly_reversible, odes, initial_conditions, chemostatted_values=None, n_samples=200, seed=0, tol=1e-09)

Attempt a Horn–Jackson global-stability certificate for one compatibility class.

odes is the species → SymPy RHS mapping (numeric rates substituted). Returns a :class:StabilityCertificate.

Source code in mantis/stability.py

def certify_global_stability(reactions, species, rate_values, deficiency,
                             is_weakly_reversible, odes, initial_conditions,
                             chemostatted_values=None, n_samples=200, seed=0,
                             tol=1e-9):
    """Attempt a Horn–Jackson global-stability certificate for one compatibility class.

    ``odes`` is the species → SymPy RHS mapping (numeric rates substituted).
    Returns a :class:`StabilityCertificate`.
    """
    import sympy
    from .analysis import _conservation_law_vectors, _make_feasible_initial

    chem = chemostatted_values or {}
    cb = is_complex_balanced(reactions, species, rate_values, deficiency,
                             is_weakly_reversible, initial_conditions, chem or None)

    if not cb:
        reason = ("not weakly reversible" if not is_weakly_reversible
                  else "not complex-balanced at these rate constants")
        lines = [
            "Global stability (Horn–Jackson Lyapunov): NOT certified",
            f"  → Network is {reason}; the complex-balanced Lyapunov function does "
            f"not apply.",
            "  → Local eigenvalue stability (see steady_states) still holds where computed.",
        ]
        return StabilityCertificate(
            is_complex_balanced=False, globally_stable=False,
            rationale=reason, summary_lines=lines,
        )

    c_star_vec = _positive_equilibrium(reactions, species, rate_values,
                                       initial_conditions, chem or None)
    if c_star_vec is None:
        lines = ["Global stability: complex-balanced, but no positive equilibrium "
                 "was found in the given class (check initial conditions)."]
        return StabilityCertificate(
            is_complex_balanced=True, globally_stable=False,
            rationale="no positive equilibrium found", summary_lines=lines,
        )
    c_star = dict(zip(species, c_star_vec.tolist()))

    V, dV = _lyapunov_symbolics(species, odes, c_star)

    # Numerically verify dV/dt ≤ 0 across the compatibility class.
    syms = [sympy.Symbol(s) for s in species]
    dV_fn = sympy.lambdify(syms, dV, "numpy")
    cl_vecs = _conservation_law_vectors(reactions, species)
    y_ref = np.array([initial_conditions.get(s, 0.0) for s in species])
    cl_totals = [float(v @ y_ref) for v in cl_vecs]
    rng = np.random.default_rng(seed)

    max_dVdt = -np.inf
    n_ok = 0
    for _ in range(n_samples):
        if cl_vecs:
            y = _make_feasible_initial(cl_vecs, cl_totals, len(species), rng)
        else:
            y = np.abs(rng.normal(loc=np.maximum(c_star_vec, 1e-9),
                                  scale=np.maximum(c_star_vec, 1e-9)))
        if y is None:
            continue
        y = np.maximum(y, 1e-30)
        val = float(dV_fn(*y))
        if not np.isfinite(val):
            continue
        max_dVdt = max(max_dVdt, val)
        n_ok += 1

    max_dVdt = 0.0 if max_dVdt == -np.inf else max_dVdt
    # Allow a small positive tolerance for floating-point noise near c*.
    numerically_ok = max_dVdt <= tol * (abs(float(np.max(c_star_vec))) + 1.0)

    lines = [
        "Global stability (Horn–Jackson Lyapunov): CERTIFIED",
        f"  → Network is complex-balanced "
        f"({'δ=0 weakly reversible, all rates' if deficiency == 0 else 'verified at equilibrium'}).",
        "  → V(c) = Σ_i [c_i(ln(c_i/c_i*) − 1) + c_i*] is a strict Lyapunov function;",
        "    dV/dt ≤ 0 with equality only at c*  ⟹  c* is GLOBALLY asymptotically",
        "    stable within its positive stoichiometric compatibility class.",
        f"  → Numerical check: dV/dt ≤ 0 at {n_ok} sampled points "
        f"(max dV/dt = {max_dVdt:.3e}).",
    ]
    return StabilityCertificate(
        is_complex_balanced=True,
        globally_stable=bool(numerically_ok),
        equilibrium=c_star,
        lyapunov_function=V,
        lyapunov_derivative=dV,
        n_samples_checked=n_ok,
        max_dVdt=max_dVdt,
        rationale="complex-balanced; Horn–Jackson Lyapunov certificate",
        summary_lines=lines,
    )

complex_balanced_equilibrium(reactions, species, rate_values, deficiency, is_weakly_reversible, initial_conditions, chemostatted_values=None)

Return the Birch point: the unique positive complex-balanced equilibrium.

For a complex-balanced network the Birch theorem guarantees exactly one positive equilibrium in each stoichiometric compatibility class (fixed here by the initial conditions). Raises ValueError if the network is not complex-balanced or no positive equilibrium is found.

Source code in mantis/stability.py

def complex_balanced_equilibrium(reactions, species, rate_values, deficiency,
                                 is_weakly_reversible, initial_conditions,
                                 chemostatted_values=None):
    """Return the Birch point: the unique positive complex-balanced equilibrium.

    For a complex-balanced network the Birch theorem guarantees exactly one positive
    equilibrium in each stoichiometric compatibility class (fixed here by the initial
    conditions).  Raises ``ValueError`` if the network is not complex-balanced or no
    positive equilibrium is found.
    """
    chem = chemostatted_values or {}
    if not is_complex_balanced(reactions, species, rate_values, deficiency,
                               is_weakly_reversible, initial_conditions, chem or None):
        raise ValueError(
            "Network is not complex-balanced at these rate constants; the Birch point "
            "is only defined for complex-balanced networks."
        )
    c_star = _positive_equilibrium(reactions, species, rate_values,
                                   initial_conditions, chem or None)
    if c_star is None:
        raise ValueError("No positive equilibrium found in the given compatibility class.")
    return dict(zip(species, c_star.tolist()))

is_complex_balanced(reactions, species, rate_values, deficiency, is_weakly_reversible, initial_conditions=None, chemostatted_values=None, tol=1e-06)

Decide whether the network is complex-balanced.

Complex balance ⇒ weak reversibility (necessary). Every weakly reversible deficiency-zero network is complex-balanced for all rates (Horn–Jackson). For weakly reversible higher-deficiency networks the property is rate-dependent and is checked numerically at the positive equilibrium (requires ICs).

Source code in mantis/stability.py

def is_complex_balanced(reactions, species, rate_values, deficiency,
                        is_weakly_reversible, initial_conditions=None,
                        chemostatted_values=None, tol=1e-6):
    """Decide whether the network is complex-balanced.

    Complex balance ⇒ weak reversibility (necessary).  Every weakly reversible
    deficiency-zero network is complex-balanced for all rates (Horn–Jackson).
    For weakly reversible higher-deficiency networks the property is rate-dependent
    and is checked numerically at the positive equilibrium (requires ICs).
    """
    if not is_weakly_reversible:
        return False
    if deficiency == 0:
        return True
    if initial_conditions is None:
        return False
    chem = chemostatted_values or {}
    c_star = _positive_equilibrium(reactions, species, rate_values,
                                   initial_conditions, chem or None)
    if c_star is None:
        return False
    imbalance = _complex_flux_balance(reactions, species, rate_values, c_star, chem, tol)
    scale = max(abs(float(rate_values.get(r.rate_key, 0.0))) for r in reactions) or 1.0
    return imbalance <= tol * scale

Absolute Concentration Robustness

mantis.acr

Absolute Concentration Robustness (ACR) detection — Shinar & Feinberg, Science 2010.

A mass-action system has absolute concentration robustness in a species S when the steady-state concentration of S is identical in every positive steady state, regardless of the initial totals (the stoichiometric compatibility class). This is exactly the robustness property a biosensor or homeostatic module is designed to have, which makes it directly relevant to diagnostic-circuit design.

Shinar & Feinberg give a purely structural sufficient condition:

If a mass-action system has deficiency one, admits a positive steady state, and
possesses two **non-terminal** complexes that differ only in species *S*, then the
system has ACR in *S*.

"Non-terminal" complexes are those that do not lie in a terminal strong linkage class (a strongly connected component of the reaction graph with no reactions leading out of it). Two complexes "differ only in S" when their species vectors differ in exactly one coordinate.

This module implements that test on top of the complex graph mantis already builds.

ACRResult dataclass

Outcome of an ACR analysis.

Attributes

has_acr : bool Whether the Shinar–Feinberg sufficient condition detected ACR. species : list[str] Species for which ACR is certified (empty when has_acr is False). deficiency : int Network deficiency (the test applies only when this is 1). applies : bool Whether the deficiency-one precondition is met. acr_values : dict[str, float] Steady-state value of each ACR species, when a positive steady state was computed (the value is independent of the initial totals). nonterminal_complex_pairs : list[tuple[str, str, str]] Witnessing (complex_a, complex_b, species) triples. summary_lines : list[str]

Source code in mantis/acr.py

@dataclass
class ACRResult:
    """Outcome of an ACR analysis.

    Attributes
    ----------
    has_acr : bool
        Whether the Shinar–Feinberg sufficient condition detected ACR.
    species : list[str]
        Species for which ACR is certified (empty when ``has_acr`` is False).
    deficiency : int
        Network deficiency (the test applies only when this is 1).
    applies : bool
        Whether the deficiency-one precondition is met.
    acr_values : dict[str, float]
        Steady-state value of each ACR species, when a positive steady state was
        computed (the value is independent of the initial totals).
    nonterminal_complex_pairs : list[tuple[str, str, str]]
        Witnessing ``(complex_a, complex_b, species)`` triples.
    summary_lines : list[str]
    """
    has_acr: bool
    species: list[str] = field(default_factory=list)
    deficiency: int = 0
    applies: bool = False
    acr_values: dict[str, float] = field(default_factory=dict)
    nonterminal_complex_pairs: list[tuple[str, str, str]] = field(default_factory=list)
    summary_lines: list[str] = field(default_factory=list)

    def __str__(self) -> str:
        return "\n".join(self.summary_lines)

detect_acr(reactions, species, graph, deficiency, rate_values=None, initial_conditions=None, chemostatted_values=None)

Detect ACR via the Shinar–Feinberg deficiency-one sufficient condition.

graph is the (reduced) complex digraph; deficiency the network deficiency. When rate_values and initial_conditions are supplied and a positive steady state exists, the robust steady-state value of each ACR species is also reported.

Source code in mantis/acr.py

def detect_acr(
    reactions: list[Reaction],
    species: list[str],
    graph: nx.DiGraph,
    deficiency: int,
    rate_values: dict[str, float] | None = None,
    initial_conditions: dict[str, float] | None = None,
    chemostatted_values: dict[str, float] | None = None,
) -> ACRResult:
    """Detect ACR via the Shinar–Feinberg deficiency-one sufficient condition.

    ``graph`` is the (reduced) complex digraph; ``deficiency`` the network deficiency.
    When ``rate_values`` and ``initial_conditions`` are supplied and a positive steady
    state exists, the robust steady-state value of each ACR species is also reported.
    """
    if deficiency != 1:
        return ACRResult(
            has_acr=False,
            deficiency=deficiency,
            applies=False,
            summary_lines=[
                "ACR (Shinar–Feinberg deficiency-one test): NOT applicable "
                f"(δ = {deficiency} ≠ 1).",
                "  → This sufficient condition only covers deficiency-one networks; "
                "absence here does not rule out ACR by other mechanisms.",
            ],
        )

    terminal = terminal_strong_linkage_classes(graph)
    terminal_nodes: set[Complex] = set().union(*terminal) if terminal else set()
    nonterminal = [c for c in graph.nodes() if c not in terminal_nodes]

    acr_species: list[str] = []
    pairs: list[tuple[str, str, str]] = []
    for i in range(len(nonterminal)):
        for j in range(i + 1, len(nonterminal)):
            ca, cb = nonterminal[i], nonterminal[j]
            diff = _complex_vector(ca, species) - _complex_vector(cb, species)
            nz = np.flatnonzero(np.abs(diff) > 1e-12)
            if len(nz) == 1:
                s = species[int(nz[0])]
                if s not in acr_species:
                    acr_species.append(s)
                pairs.append((_complex_label(ca), _complex_label(cb), s))

    has_acr = len(acr_species) > 0

    acr_values: dict[str, float] = {}
    if has_acr and rate_values is not None and initial_conditions is not None:
        from .analysis import find_steady_states
        states = find_steady_states(
            reactions, species, rate_values, initial_conditions,
            n_attempts=12, chemostatted_values=chemostatted_values or None,
        )
        for st in states:
            y = np.array([st.concentrations[s] for s in species])
            if np.all(y >= 0) and st.residual < 1e-4 * (np.max(np.abs(y)) + 1e-30):
                for s in acr_species:
                    acr_values.setdefault(s, float(st.concentrations[s]))
                break

    if has_acr:
        lines = [
            "ACR (Shinar–Feinberg deficiency-one test): DETECTED",
            f"  → Robust species: {', '.join(sorted(acr_species))}",
            "  → Their steady-state concentrations are independent of the initial "
            "totals (within any class admitting a positive steady state).",
        ]
        for a, b, s in pairs:
            lines.append(f"     witness: non-terminal complexes ‘{a}’ and ‘{b}’ differ only in {s}")
        for s, v in sorted(acr_values.items()):
            lines.append(f"     {s}* = {v:.6g} (computed at a positive steady state)")
    else:
        lines = [
            "ACR (Shinar–Feinberg deficiency-one test): NOT detected",
            "  → δ = 1 but no two non-terminal complexes differ in a single species.",
        ]

    return ACRResult(
        has_acr=has_acr,
        species=sorted(acr_species),
        deficiency=deficiency,
        applies=True,
        acr_values=acr_values,
        nonterminal_complex_pairs=pairs,
        summary_lines=lines,
    )

terminal_strong_linkage_classes(graph)

Return the terminal strong linkage classes of a complex graph.

A strong linkage class (SCC) is terminal when no edge leaves it for another SCC — i.e. it is a sink in the condensation.

Source code in mantis/acr.py

def terminal_strong_linkage_classes(graph: nx.DiGraph) -> list[set[Complex]]:
    """Return the terminal strong linkage classes of a complex graph.

    A strong linkage class (SCC) is *terminal* when no edge leaves it for another
    SCC — i.e. it is a sink in the condensation.
    """
    sccs = list(nx.strongly_connected_components(graph))
    condensation = nx.condensation(graph, scc=sccs)
    terminal = []
    for node in condensation.nodes():
        if condensation.out_degree(node) == 0:
            terminal.append(set(condensation.nodes[node]["members"]))
    return terminal

Injectivity

mantis.injectivity

Injectivity test — Craciun & Feinberg (SIAM J. Appl. Math. 2005, 2006).

A mass-action network is injective when its species-formation-rate map is injective on each stoichiometric compatibility class; an injective network has at most one positive steady state per class, so it cannot be multistationary. Crucially this rules out multistationarity for higher-deficiency networks too, where the Deficiency One Theorem is silent.

Craciun & Feinberg reduce injectivity to a sign condition on a determinant. Form the square "steady-state map"

F(c) = ( rank(N) independent components of N·v(c) ,  the conservation totals ),

whose Jacobian J(c) is an n×n matrix whose entries are polynomials in the (positive) concentrations c and the (positive) rate constants k. If det J(c) is sign-definite — every monomial coefficient shares one sign, so the determinant never vanishes for any positive c and k — then F is injective and the network admits at most one positive steady state per class.

The test is sufficient: a sign-definite determinant certifies injectivity, while a sign-indefinite one is inconclusive (the network may still be injective for specific rates, or may genuinely be multistationary, as for the Schlögl model).

InjectivityResult dataclass

Outcome of the Craciun–Feinberg injectivity test.

Attributes

injective : bool True ⇒ the network is certified injective (at most one positive steady state per compatibility class). False ⇒ inconclusive (not sign-definite). is_sign_definite : bool Whether the steady-state Jacobian determinant is sign-definite. determinant : sympy.Expr The symbolic determinant (in concentrations and rate constants). reason : str summary_lines : list[str]

Source code in mantis/injectivity.py

@dataclass
class InjectivityResult:
    """Outcome of the Craciun–Feinberg injectivity test.

    Attributes
    ----------
    injective : bool
        True ⇒ the network is certified injective (at most one positive steady state
        per compatibility class). False ⇒ inconclusive (not sign-definite).
    is_sign_definite : bool
        Whether the steady-state Jacobian determinant is sign-definite.
    determinant : sympy.Expr
        The symbolic determinant (in concentrations and rate constants).
    reason : str
    summary_lines : list[str]
    """
    injective: bool
    is_sign_definite: bool
    determinant: "object" = None
    reason: str = ""
    summary_lines: list[str] = field(default_factory=list)

    def __str__(self) -> str:
        return "\n".join(self.summary_lines)

test_injectivity(reactions, species, chemostatted_values=None)

Run the determinant sign-definiteness injectivity test (symbolic, rate-free).

Source code in mantis/injectivity.py

def test_injectivity(
    reactions: list[Reaction],
    species: list[str],
    chemostatted_values: dict[str, float] | None = None,
) -> InjectivityResult:
    """Run the determinant sign-definiteness injectivity test (symbolic, rate-free)."""
    import sympy

    from .symbolic import make_species_symbols, make_rate_symbols, mass_action_flux

    chem = set((chemostatted_values or {}).keys())

    # Symbols for every species appearing in the reactions (dynamic + chemostatted),
    # so chemostatted species enter the determinant as positive parameters.
    all_names: list[str] = []
    seen: set[str] = set()
    for rxn in reactions:
        for name, _ in list(rxn.reactants) + list(rxn.products):
            if name not in seen:
                seen.add(name)
                all_names.append(name)
    sp_syms = make_species_symbols(sorted(all_names))
    rate_syms, key_to_sym = make_rate_symbols(reactions)

    dyn_syms = [sp_syms[s] for s in species]
    n = len(species)
    if n == 0:
        return InjectivityResult(True, True, sympy.Integer(0), "no dynamic species",
                                 ["Injectivity: trivially injective (no dynamic species)."])

    # Symbolic mass-action fluxes (all species symbolic — nothing folded to a number).
    fluxes = [mass_action_flux(rxn, sp_syms, key_to_sym[rxn.rate_key], None)
              for rxn in reactions]

    N = build_stoichiometry_matrix(reactions, species)
    rank = matrix_rank(N)

    # rank(N) independent rows of N·v, chosen by QR column pivoting on Nᵀ.
    if rank > 0:
        from scipy.linalg import qr
        _, _, perm = qr(N.T, pivoting=True)
        ind = list(perm[:rank])
    else:
        ind = []

    rows = []
    for i in ind:
        expr = sympy.Integer(0)
        for j in range(len(reactions)):
            coeff = int(round(N[i, j]))
            if coeff:
                expr += coeff * fluxes[j]
        rows.append(sympy.expand(expr))

    # Conservation-law rows (their gradients pin down the remaining n - rank directions).
    laws = conservation_laws_sympy(N, species)
    for law in laws:
        rows.append(law)

    if len(rows) != n:
        return InjectivityResult(
            False, False, None,
            "could not form a square steady-state map",
            ["Injectivity: inconclusive (degenerate stoichiometry)."],
        )

    J = sympy.Matrix([[sympy.diff(r, c) for c in dyn_syms] for r in rows])
    det = sympy.expand(J.det())

    if det == 0:
        return InjectivityResult(
            False, False, det,
            "determinant is identically zero",
            [
                "Injectivity (Craciun–Feinberg): NOT certified",
                "  → The steady-state Jacobian determinant vanishes identically "
                "(degenerate); the test is inconclusive.",
            ],
        )

    poly = sympy.Poly(det, *(dyn_syms + list(rate_syms.values())
                             + [sp_syms[s] for s in sorted(all_names) if s in chem]))
    signs = {int(np.sign(float(c))) for c in poly.coeffs() if c != 0}
    sign_definite = len(signs) == 1

    if sign_definite:
        lines = [
            "Injectivity (Craciun–Feinberg): CERTIFIED injective",
            "  → The steady-state Jacobian determinant is sign-definite, so it never "
            "vanishes for positive concentrations and rate constants.",
            "  → The network admits AT MOST ONE positive steady state per stoichiometric",
            "    compatibility class — multistationarity is structurally impossible",
            "    (for all rate constants).",
        ]
    else:
        lines = [
            "Injectivity (Craciun–Feinberg): NOT certified (inconclusive)",
            "  → The determinant is sign-indefinite; this sufficient test cannot rule out",
            "    multiple steady states (the network may be multistationary for some rates).",
        ]

    return InjectivityResult(
        injective=sign_definite,
        is_sign_definite=sign_definite,
        determinant=det,
        reason="sign-definite determinant" if sign_definite else "sign-indefinite determinant",
        summary_lines=lines,
    )

Pseudo-arclength continuation (fold / Hopf)

mantis.continuation

Pseudo-arclength continuation of steady-state branches (fold / Hopf detection).

The :meth:~mantis.network.CRNetwork.bifurcation log-scan recomputes steady states independently at each parameter value. Near a saddle-node (fold) the branch turns back on itself, so two states share one parameter value and one disappears as the parameter is swept past the fold — a naive scan loses the branch, mis-pairs points across the fold, and never reports the fold itself.

Pseudo-arclength continuation (Keller 1977; Allgower–Georg 2003) fixes this by parameterising the branch not by the bifurcation parameter λ but by arclength s along the solution curve (c(s), λ(s)). Folds — where dλ/ds = 0 — are then ordinary points of the curve and are traversed smoothly. The method is the standard predictor–corrector:

1. Predictor — step ds along the unit tangent of the branch.

2. Corrector — Newton on the augmented square system

   G(c, λ) = 0 (n steady-state equations) ⟨t_c, c − c_prev⟩ + t_λ (λ − λ_prev) − ds = 0 (arclength constraint)

The arclength row keeps the Jacobian non-singular through the fold, where the steady-state Jacobian G_c alone is singular.

G is the conservation-constrained steady-state map: rank(N) independent rows of the mass-action ODE field stacked on the n − rank(N) conservation laws (totals fixed by the initial conditions), so the branch stays inside one stoichiometric compatibility class.

Bifurcation detection along the branch:

• Fold (saddle-node) — det G_c changes sign (the steady-state Jacobian is singular exactly at a fold), the canonical fold test function.
• Hopf — the number of eigenvalues with positive real part changes without a fold (so a complex-conjugate pair, not a real eigenvalue, has crossed the imaginary axis). Reported as a candidate; a Hopf marks the birth of a limit cycle. Eigenvalues are those of the Jacobian restricted to the stoichiometric subspace (the n − rank(N) exact zeros from conservation are dropped).

The classic demonstration is the Schlögl model: continuation traces the full S-shaped branch and pins both fold points bounding the bistable region, neither of which the log-scan reports.

BifurcationPoint dataclass

A detected bifurcation on a continuation branch.

Attributes

kind : str "fold" (saddle-node) or "hopf" (candidate Hopf). parameter : float Bifurcation-parameter value at the point (linearly localised between the two bracketing continuation steps). state : dict[str, float] Species concentrations at the bifurcation. eigenvalue : complex The critical eigenvalue nearest the imaginary axis at the crossing.

Source code in mantis/continuation.py

@dataclass
class BifurcationPoint:
    """A detected bifurcation on a continuation branch.

    Attributes
    ----------
    kind : str
        ``"fold"`` (saddle-node) or ``"hopf"`` (candidate Hopf).
    parameter : float
        Bifurcation-parameter value at the point (linearly localised between the
        two bracketing continuation steps).
    state : dict[str, float]
        Species concentrations at the bifurcation.
    eigenvalue : complex
        The critical eigenvalue nearest the imaginary axis at the crossing.
    """
    kind: str
    parameter: float
    state: dict[str, float]
    eigenvalue: complex

    def __str__(self) -> str:
        loc = ", ".join(f"{s}={v:.4g}" for s, v in self.state.items())
        return (f"{self.kind.upper()} at {self.parameter:.6g} "
                f"(λ_crit={self.eigenvalue:.3g}); {loc}")

ContinuationResult dataclass

A steady-state branch traced by pseudo-arclength continuation.

The branch is ordered by arclength, so parameter_values is generally non-monotonic (it reverses at each fold) — that is the whole point.

Attributes

parameter_name : str Rate key that was continued. parameter_values : np.ndarray Bifurcation-parameter value at each branch point (arclength order). branch : list[dict[str, float]] Species concentrations at each branch point. stable : np.ndarray Boolean stability of each branch point (all reduced eigenvalues Re < 0). eigenvalues : list[np.ndarray] Reduced-Jacobian eigenvalues at each branch point. bifurcations : list[BifurcationPoint]

Source code in mantis/continuation.py

@dataclass
class ContinuationResult:
    """A steady-state branch traced by pseudo-arclength continuation.

    The branch is ordered by arclength, so ``parameter_values`` is generally
    *non-monotonic* (it reverses at each fold) — that is the whole point.

    Attributes
    ----------
    parameter_name : str
        Rate key that was continued.
    parameter_values : np.ndarray
        Bifurcation-parameter value at each branch point (arclength order).
    branch : list[dict[str, float]]
        Species concentrations at each branch point.
    stable : np.ndarray
        Boolean stability of each branch point (all reduced eigenvalues Re < 0).
    eigenvalues : list[np.ndarray]
        Reduced-Jacobian eigenvalues at each branch point.
    bifurcations : list[BifurcationPoint]
    """
    parameter_name: str
    parameter_values: np.ndarray
    branch: list[dict[str, float]]
    stable: np.ndarray
    eigenvalues: list[np.ndarray] = field(default_factory=list)
    bifurcations: list[BifurcationPoint] = field(default_factory=list)

    def folds(self) -> list[BifurcationPoint]:
        return [b for b in self.bifurcations if b.kind == "fold"]

    def hopfs(self) -> list[BifurcationPoint]:
        return [b for b in self.bifurcations if b.kind == "hopf"]

    def species_branch(self, species: str) -> np.ndarray:
        """The concentration of one species along the branch (arclength order)."""
        return np.array([pt[species] for pt in self.branch])

    def __str__(self) -> str:
        lines = [
            f"Continuation in {self.parameter_name!r}: {len(self.branch)} branch "
            f"points over λ ∈ [{self.parameter_values.min():.4g}, "
            f"{self.parameter_values.max():.4g}]",
        ]
        if self.bifurcations:
            lines.append(f"  {len(self.bifurcations)} bifurcation(s):")
            lines += [f"    • {b}" for b in self.bifurcations]
        else:
            lines.append("  no bifurcations detected on this branch")
        return "\n".join(lines)

species_branch(species)

The concentration of one species along the branch (arclength order).

Source code in mantis/continuation.py

def species_branch(self, species: str) -> np.ndarray:
    """The concentration of one species along the branch (arclength order)."""
    return np.array([pt[species] for pt in self.branch])

pseudo_arclength_continuation(reactions, species, base_rates, parameter, param_range, initial_conditions, chemostatted_values=None, ds=0.05, max_steps=2000, newton_tol=1e-10, initial_state=None)

Trace a steady-state branch by pseudo-arclength continuation.

Parameters

parameter : str Rate key to continue (e.g. "A + 2X -> 3X"). param_range : (float, float) (λ_min, λ_max) box the branch is traced within; continuation stops when it leaves the box. To capture an S-curve the box must bracket both folds. initial_conditions : dict[str, float] Fixes the conservation totals (the compatibility class) and, unless initial_state is given, seeds the branch (a steady state is solved at λ_min to start from). ds : float Arclength step (in the non-dimensional (c/scale, λ/λ_min) metric); adapted down on Newton failure and back up on easy steps. max_steps : int Cap on continuation steps (guards against closed isolas / cycling). initial_state : dict[str, float], optional Explicit starting steady state on the desired branch; overrides the automatic seed (use to select which branch to follow).

Returns

ContinuationResult

Source code in mantis/continuation.py

def pseudo_arclength_continuation(
    reactions: list[Reaction],
    species: list[str],
    base_rates: dict[str, float],
    parameter: str,
    param_range: tuple[float, float],
    initial_conditions: dict[str, float],
    chemostatted_values: dict[str, float] | None = None,
    ds: float = 0.05,
    max_steps: int = 2000,
    newton_tol: float = 1e-10,
    initial_state: dict[str, float] | None = None,
) -> ContinuationResult:
    """Trace a steady-state branch by pseudo-arclength continuation.

    Parameters
    ----------
    parameter : str
        Rate key to continue (e.g. ``"A + 2X -> 3X"``).
    param_range : (float, float)
        ``(λ_min, λ_max)`` box the branch is traced within; continuation stops
        when it leaves the box.  To capture an S-curve the box must bracket both
        folds.
    initial_conditions : dict[str, float]
        Fixes the conservation totals (the compatibility class) and, unless
        ``initial_state`` is given, seeds the branch (a steady state is solved at
        ``λ_min`` to start from).
    ds : float
        Arclength step (in the non-dimensional ``(c/scale, λ/λ_min)`` metric);
        adapted down on Newton failure and back up on easy steps.
    max_steps : int
        Cap on continuation steps (guards against closed isolas / cycling).
    initial_state : dict[str, float], optional
        Explicit starting steady state on the desired branch; overrides the
        automatic seed (use to select which branch to follow).

    Returns
    -------
    ContinuationResult
    """
    from .parsing import normalize_rate_key
    from .analysis import find_steady_states, classify_steady_state
    from .steady_states import all_steady_states

    try:
        param_key = normalize_rate_key(parameter)
    except ValueError:
        param_key = parameter

    lam_min, lam_max = float(param_range[0]), float(param_range[1])
    chem_vals = chemostatted_values or {}
    smap = _SteadyStateMap(reactions, species, base_rates, param_key,
                           chem_vals, initial_conditions)
    n = smap.n

    # ── Seed: a steady state on the branch at λ_min ──────────────────────────
    if initial_state is not None:
        c = np.array([initial_state.get(s, 0.0) for s in species], float)
    else:
        seed_rates = dict(base_rates)
        seed_rates[param_key] = lam_min
        # Exhaustive enumeration gives a reliable seed (the multi-start solver can
        # converge to spurious points on open systems); fall back to it otherwise.
        ss = []
        try:
            ss = all_steady_states(reactions, species, seed_rates, initial_conditions,
                                   chemostatted_values=chem_vals or None)
        except Exception:
            ss = []
        if not ss:
            ss = find_steady_states(reactions, species, seed_rates, initial_conditions,
                                    n_attempts=20, seed=0,
                                    chemostatted_values=chem_vals or None)
        if not ss:
            raise RuntimeError("Could not find a steady state to seed continuation.")
        c = np.array([ss[0].concentrations[s] for s in species], float)
    lam = lam_min

    # Non-dimensionalising scales so c and λ contribute comparably to arclength.
    cscale = max(float(np.max(np.abs(c))), 1.0)
    lscale = max(abs(lam_min), abs(lam_max), 1e-300)

    def scaled(cv, lv):
        return np.concatenate([cv / cscale, [lv / lscale]])

    def unscaled(x):
        return x[:n] * cscale, x[n] * lscale

    # ── Record helpers ───────────────────────────────────────────────────────
    params: list[float] = []
    branch: list[dict[str, float]] = []
    stable: list[bool] = []
    eig_list: list[np.ndarray] = []
    det_signs: list[float] = []
    n_unstable: list[int] = []
    bifs: list[BifurcationPoint] = []

    xs_scaled: list[np.ndarray] = []      # accepted points in scaled coordinates
    tans: list[np.ndarray] = []           # unit tangent at each accepted point

    def record(cv, lv, tan):
        eigs = smap.reduced_eigenvalues(cv, lv)
        is_stable, _ = classify_steady_state(eigs)
        params.append(lv)
        branch.append(dict(zip(species, cv.tolist())))
        stable.append(is_stable)
        eig_list.append(eigs)
        det = float(np.linalg.det(smap.G_c(cv, lv)))
        det_signs.append(np.sign(det) if det != 0 else (det_signs[-1] if det_signs else 1.0))
        n_unstable.append(int(np.sum(eigs.real > 1e-9 * (np.max(np.abs(eigs)) + 1e-300))))
        xs_scaled.append(scaled(cv, lv))
        tans.append(tan)

    record(c, lam, None)
    prev_tan = None

    for _ in range(max_steps):
        Gc = smap.G_c(c, lam)
        Glam = smap.G_lam(c, lam)
        # Scale derivatives consistently with the (c/cscale, λ/lscale) variables.
        tan = _tangent(Gc * cscale, Glam * lscale, prev_tan)

        x0 = scaled(c, lam)
        x_pred = x0 + ds * tan

        x, converged = _correct(smap, n, cscale, lscale, x_pred, tan, x0, ds, newton_tol)

        if not converged:
            ds *= 0.5
            if ds < 1e-9:
                break
            continue

        cv, lv = unscaled(x)
        if np.any(cv < -1e-7 * (np.max(np.abs(cv)) + 1e-30)):
            break                                   # left the positive orthant
        cv = np.maximum(cv, 0.0)

        if tans[-1] is None:
            tans[-1] = tan                          # tangent leaving the seed point
        record(cv, lv, tan)

        c, lam = cv, lv
        prev_tan = tan
        ds = min(ds * 1.1, 0.1)                      # grow step on success

        if lam < lam_min or lam > lam_max:
            break

    _detect_bifurcations(smap, species, n, cscale, lscale, xs_scaled, tans,
                         params, det_signs, n_unstable, bifs)

    return ContinuationResult(
        parameter_name=param_key,
        parameter_values=np.array(params),
        branch=branch,
        stable=np.array(stable),
        eigenvalues=eig_list,
        bifurcations=bifs,
    )

Multistationarity parameter regions (Conradi–Feliu–Mincheva–Wiuf)

mantis.multistationarity

Parameter regions for multistationarity — Conradi, Feliu, Mincheva & Wiuf (PLoS Comput. Biol. 2017).

Injectivity (:mod:mantis.injectivity) answers a yes/no question — can this network ever have more than one positive steady state in a compatibility class? Conradi et al. (2017) go further and ask where in parameter space it does: they give the determinant criterion that carves the rate-/total-concentration space into a monostationary region and a multistationary region.

The object of study is the critical function

φ(c, k) = (−1)^s · det J(c, k),

where J is the Jacobian of the canonical steady-state map (s = rank(N) independent mass-action ODE rows stacked on the n − s conservation laws) and s is the dimension of the stoichiometric subspace. Reading φ as a polynomial in the steady-state concentrations c with coefficients that are polynomials in the rate constants k:

• If every coefficient is positive for all positive k then φ > 0 on the whole positive orthant, the steady-state map is injective, and the network is monostationary for all rates and totals (Craciun–Feinberg).
• If some coefficient can be negative, φ attains negative values and the network is multistationary for parameters in that region. Each such coefficient inequality a_α(k) < 0 is a face of the multistationarity region — the explicit "where" that the 2017 paper contributes.

For a network whose positive steady states admit a monomial parametrisation (weakly reversible / toric steady states, the paper's setting) a negative coefficient is both necessary and sufficient for multistationarity; in general it is the necessary condition (a sign-change witness). This module returns the critical function, the per-coefficient sign certificate, and — when numeric rates are supplied — whether those specific rates land in the multistationary region.

MultistationarityResult dataclass

Outcome of the Conradi–Feliu–Mincheva–Wiuf critical-function analysis.

Attributes

monostationary : bool True ⇒ certified at most one positive steady state per compatibility class for all positive rate constants (the critical function is sign-definite positive). multistationary_possible : bool True ⇒ some coefficient of the critical function can be negative, so multistationarity occurs for parameters in region_conditions. critical_function : sympy.Expr φ = (−1)^s · det J in concentrations and rate constants. coefficient_terms : list[tuple] (monomial, coefficient_expr, sign) for each concentration-monomial of φ; sign ∈ {"+", "−", "±"} is the coefficient's sign over positive rate constants. region_conditions : list[sympy.Expr] Coefficient expressions whose negativity delimits the multistationarity region (expr < 0 ⇒ multistationarity). Empty when monostationary. multistationary_at_rates : bool | None If numeric rates were supplied: whether they actually yield ≥ 2 positive steady states. When the steady-state set reduces to a univariate polynomial this is exact (a positive-root count); otherwise it falls back to the necessary sign condition (a critical coefficient is negative). None if no rates were supplied. region_boundary : sympy.Expr | None The exact multistationarity-region boundary when the steady states reduce to a univariate polynomial p(x): the discriminant disc_x(p). Crossing region_boundary = 0 changes the number of (positive) steady states, so it is the genuine fold locus in parameter space — unlike region_conditions, which is only the necessary sign condition. None when the reduction is unavailable. steady_state_polynomial : sympy.Expr | None The reduced univariate steady-state polynomial p(x) (in one species, with rate / total parameters), when available. summary_lines : list[str]

Source code in mantis/multistationarity.py

@dataclass
class MultistationarityResult:
    """Outcome of the Conradi–Feliu–Mincheva–Wiuf critical-function analysis.

    Attributes
    ----------
    monostationary : bool
        True ⇒ certified at most one positive steady state per compatibility
        class for *all* positive rate constants (the critical function is
        sign-definite positive).
    multistationary_possible : bool
        True ⇒ some coefficient of the critical function can be negative, so
        multistationarity occurs for parameters in ``region_conditions``.
    critical_function : sympy.Expr
        ``φ = (−1)^s · det J`` in concentrations and rate constants.
    coefficient_terms : list[tuple]
        ``(monomial, coefficient_expr, sign)`` for each concentration-monomial of
        ``φ``; ``sign`` ∈ {``"+"``, ``"−"``, ``"±"``} is the coefficient's sign
        over positive rate constants.
    region_conditions : list[sympy.Expr]
        Coefficient expressions whose negativity delimits the multistationarity
        region (``expr < 0`` ⇒ multistationarity).  Empty when monostationary.
    multistationary_at_rates : bool | None
        If numeric rates were supplied: whether they actually yield ≥ 2 positive
        steady states.  When the steady-state set reduces to a univariate
        polynomial this is *exact* (a positive-root count); otherwise it falls
        back to the necessary sign condition (a critical coefficient is negative).
        ``None`` if no rates were supplied.
    region_boundary : sympy.Expr | None
        The **exact** multistationarity-region boundary when the steady states
        reduce to a univariate polynomial ``p(x)``: the discriminant
        ``disc_x(p)``.  Crossing ``region_boundary = 0`` changes the number of
        (positive) steady states, so it is the genuine fold locus in parameter
        space — unlike ``region_conditions``, which is only the necessary sign
        condition.  ``None`` when the reduction is unavailable.
    steady_state_polynomial : sympy.Expr | None
        The reduced univariate steady-state polynomial ``p(x)`` (in one species,
        with rate / total parameters), when available.
    summary_lines : list[str]
    """
    monostationary: bool
    multistationary_possible: bool
    critical_function: object = None
    coefficient_terms: list = field(default_factory=list)
    region_conditions: list = field(default_factory=list)
    multistationary_at_rates: "bool | None" = None
    summary_lines: list[str] = field(default_factory=list)
    region_boundary: object = None
    steady_state_polynomial: object = None

    def __str__(self) -> str:
        return "\n".join(self.summary_lines)

multistationarity_region(reactions, species, rate_values=None, chemostatted_values=None)

Compute the multistationarity parameter region via the CFMW critical function.

Source code in mantis/multistationarity.py

def multistationarity_region(
    reactions: list[Reaction],
    species: list[str],
    rate_values: dict[str, float] | None = None,
    chemostatted_values: dict[str, float] | None = None,
) -> MultistationarityResult:
    """Compute the multistationarity parameter region via the CFMW critical function."""
    import sympy

    from .symbolic import make_species_symbols, make_rate_symbols, mass_action_flux

    chem = chemostatted_values or {}
    chem_keys = set(chem)

    # Symbols for every species in the reactions; chemostatted ones enter φ as
    # positive parameters (never folded to a number, so they appear in the region).
    all_names: list[str] = []
    seen: set[str] = set()
    for rxn in reactions:
        for name, _ in list(rxn.reactants) + list(rxn.products):
            if name not in seen:
                seen.add(name)
                all_names.append(name)
    sp_syms = make_species_symbols(sorted(all_names))
    rate_syms, key_to_sym = make_rate_symbols(reactions)

    dyn_syms = [sp_syms[s] for s in species]
    n = len(species)
    N = build_stoichiometry_matrix(reactions, species)
    s = matrix_rank(N)

    if n == 0:
        return MultistationarityResult(
            True, False, sympy.Integer(0), [], [], None,
            ["Multistationarity (Conradi–Feliu–Mincheva–Wiuf): no dynamic species."],
        )

    fluxes = [mass_action_flux(rxn, sp_syms, key_to_sym[rxn.rate_key], None)
              for rxn in reactions]

    # rank(N) independent ODE rows (QR column pivot on Nᵀ) + the conservation laws.
    if s > 0:
        from scipy.linalg import qr
        _, _, perm = qr(N.T, pivoting=True)
        ind = list(perm[:s])
    else:
        ind = []

    rows = []
    for i in ind:
        expr = sympy.Integer(0)
        for j in range(len(reactions)):
            coeff = int(round(N[i, j]))
            if coeff:
                expr += coeff * fluxes[j]
        rows.append(sympy.expand(expr))
    rows += list(conservation_laws_sympy(N, species))

    if len(rows) != n:
        return MultistationarityResult(
            False, True, None, [], [], None,
            ["Multistationarity: inconclusive (degenerate stoichiometry)."],
        )

    J = sympy.Matrix([[sympy.diff(r, c) for c in dyn_syms] for r in rows])
    det = sympy.expand(J.det())
    phi = sympy.expand((-1) ** s * det)

    # Positive parameters: rate constants and chemostatted concentrations.
    chem_param_syms = [sp_syms[s_] for s_ in sorted(all_names) if s_ in chem_keys]
    positive_syms = list(rate_syms.values()) + chem_param_syms

    if phi == 0:
        return MultistationarityResult(
            False, True, phi, [], [], None,
            [
                "Multistationarity (Conradi–Feliu–Mincheva–Wiuf): INCONCLUSIVE",
                "  → The critical function vanishes identically (degenerate Jacobian).",
            ],
        )

    # φ as a polynomial in the concentrations; group coefficients (in k, chemostats).
    conc_poly = sympy.Poly(phi, *dyn_syms) if any(c in phi.free_symbols for c in dyn_syms) \
        else None
    coeff_terms = []
    region = []
    if conc_poly is not None:
        for monom, coeff in conc_poly.terms():
            mon_expr = sympy.prod([c ** e for c, e in zip(dyn_syms, monom)])
            sign = _coefficient_sign(coeff, positive_syms)
            coeff_terms.append((mon_expr, sympy.expand(coeff), sign))
            if sign in ("−", "±"):
                region.append(sympy.expand(coeff))
    else:
        sign = _coefficient_sign(phi, positive_syms)
        coeff_terms.append((sympy.Integer(1), phi, sign))
        if sign in ("−", "±"):
            region.append(phi)

    all_positive = all(t[2] in ("+", "0") for t in coeff_terms)
    monostationary = all_positive
    multistationary_possible = not all_positive

    # ── Exact region: restrict to the steady-state variety (univariate reduction) ─
    # Reading φ's coefficients over the whole positive orthant is only *necessary*
    # for multistationarity (φ < 0 may occur off the steady-state set).  When the
    # steady states reduce to a univariate polynomial p(x), multistationarity ⇔ p
    # has ≥ 2 positive roots, whose boundary is the discriminant disc_x(p) = 0 — the
    # exact fold locus.
    region_boundary = None
    ss_poly = None
    var = None
    try:
        poly_info = _steady_state_polynomial(
            reactions, species, sp_syms, rate_syms, key_to_sym, dyn_syms, N, s
        )
    except Exception:
        poly_info = None
    if poly_info is not None:
        p, var, _totals = poly_info
        if p.degree() >= 2:
            ss_poly = p.as_expr()
            try:
                region_boundary = sympy.factor(sympy.discriminant(p, var))
            except Exception:
                region_boundary = None

    # ── Locate the supplied numeric rates relative to the region ──────────────
    at_rates = None
    rate_subs = None
    if rate_values is not None:
        rate_subs = {key_to_sym[k]: sympy.Float(v)
                     for k, v in rate_values.items() if k in key_to_sym}
        rate_subs.update({sp_syms[s_]: sympy.Float(chem[s_])
                          for s_ in chem_keys if s_ in sp_syms})

    # Exact numeric verdict via positive-root count, when the polynomial is fully
    # determined by the supplied rates (no free conservation totals remain).
    exact_at_rates = False
    if rate_subs is not None and ss_poly is not None:
        npos = _count_positive_roots(p, var, rate_subs)
        if npos is not None:
            at_rates = npos >= 2
            exact_at_rates = True

    # Fallback: the necessary sign condition (a critical coefficient goes negative).
    if not exact_at_rates and rate_subs is not None and multistationary_possible:
        at_rates = False
        for _mon, coeff, _sign in coeff_terms:
            val = coeff.subs(rate_subs)
            if val.free_symbols:
                continue
            if float(val) < 0:
                at_rates = True
                break

    # ── Human-readable summary ───────────────────────────────────────────────
    if monostationary:
        lines = [
            "Multistationarity (Conradi–Feliu–Mincheva–Wiuf): MONOSTATIONARY",
            f"  → The critical function φ = (−1)^{s}·det J is sign-definite positive,",
            "    so it never vanishes for positive concentrations and rate constants.",
            "  → At most ONE positive steady state per compatibility class, for ALL",
            "    rate constants and total concentrations — multistationarity is",
            "    structurally impossible.",
        ]
    else:
        lines = [
            "Multistationarity (Conradi–Feliu–Mincheva–Wiuf): REGION EXISTS",
            f"  → The critical function φ = (−1)^{s}·det J has coefficient(s) that can",
            "    be negative, so multistationarity occurs in a region of parameter space.",
        ]
        if region_boundary is not None:
            # Exact boundary on the steady-state variety (univariate reduction).
            lines.append("  → Exact region boundary (steady states reduce to a univariate")
            lines.append(f"    polynomial in {var}): the discriminant")
            lines.append(f"        disc({var}) = {region_boundary}")
            lines.append("    vanishes at the folds; the sign of disc fixes the number of")
            lines.append("    (positive) steady states on either side.")
        else:
            lines.append("  → Necessary sign condition (φ attains a negative value) requires:")
            for coeff in region[:8]:
                lines.append(f"        {sympy.simplify(coeff)} < 0")
            if len(region) > 8:
                lines.append(f"        … and {len(region) - 8} more coefficient condition(s)")
        if at_rates is not None:
            if exact_at_rates:
                verdict = ("multistationary (≥ 2 positive steady states)"
                           if at_rates else
                           "monostationary (one positive steady state)")
                lines.append(f"  → At the supplied rates the network is {verdict} "
                             "— exact positive-root count.")
            else:
                verdict = ("satisfy the necessary sign condition (multistationarity "
                           "occurs for some total concentrations)" if at_rates else
                           "do NOT make any critical coefficient negative on their own")
                lines.append(f"  → The supplied rate constants {verdict}.")

    return MultistationarityResult(
        monostationary=monostationary,
        multistationary_possible=multistationary_possible,
        critical_function=phi,
        coefficient_terms=coeff_terms,
        region_conditions=region,
        multistationary_at_rates=at_rates,
        summary_lines=lines,
        region_boundary=region_boundary,
        steady_state_polynomial=ss_poly,
    )

Stochastic stationary distribution (Anderson–Craciun–Kurtz)

mantis.stochastic_stationary

Exact stochastic stationary distributions — Anderson, Craciun & Kurtz (Bull. Math. Biol. 2010).

For a complex-balanced mass-action network the stochastic (chemical master equation) model has an exact stationary distribution in closed form — a product of independent Poissons whose means are the deterministic complex-balanced equilibrium scaled to molecule counts:

π(x) = M · Π_i  λ_i^{x_i} / x_i!,        λ_i = (N_A · V) · c_i*,

supported on the reachable state space (the stoichiometric compatibility class). With no conservation laws the state space is all of ℤ≥0ⁿ and π factorises into independent Poisson(λ_i); with conservation laws π is that product conditioned on the conserved totals (for a single moiety this is exactly a Binomial / multinomial). This lets you write the stationary distribution down from network structure alone — no simulation — for any weakly reversible deficiency-zero network, the regime where SSA is noisiest.

StationaryDistribution dataclass

Exact product-form stationary distribution of a complex-balanced network.

Attributes

species : list[str] means : np.ndarray The Poisson means λ_i (in molecule counts) = N_A · V · c_i*. volume_L : float is_conditioned : bool Whether conservation laws restrict the support to a compatibility class. is_exact : bool True when the distribution is represented exactly (independent Poissons, or a fully enumerated finite class); False when the class could not be enumerated and the independent-Poisson factorisation is used as an approximation.

Source code in mantis/stochastic_stationary.py

@dataclass
class StationaryDistribution:
    """Exact product-form stationary distribution of a complex-balanced network.

    Attributes
    ----------
    species : list[str]
    means : np.ndarray
        The Poisson means ``λ_i`` (in molecule counts) = ``N_A · V · c_i*``.
    volume_L : float
    is_conditioned : bool
        Whether conservation laws restrict the support to a compatibility class.
    is_exact : bool
        True when the distribution is represented exactly (independent Poissons, or a
        fully enumerated finite class); False when the class could not be enumerated
        and the independent-Poisson factorisation is used as an approximation.
    """
    species: list[str]
    means: np.ndarray
    volume_L: float
    is_conditioned: bool
    is_exact: bool
    _states: np.ndarray | None = field(default=None, repr=False)
    _logZ: float | None = field(default=None, repr=False)

    # ── core ────────────────────────────────────────────────────────────────────
    def poisson_means(self) -> dict[str, float]:
        """Poisson means ``λ_i`` keyed by species (molecule counts)."""
        return {s: float(m) for s, m in zip(self.species, self.means)}

    def _log_weight(self, x: np.ndarray) -> float:
        from scipy.special import gammaln
        x = np.asarray(x, dtype=float)
        lam = np.maximum(self.means, 1e-300)
        return float(np.sum(x * np.log(lam) - gammaln(x + 1.0)))

    def probability(self, counts: dict[str, float]) -> float:
        """Normalised stationary probability of an integer state ``counts``."""
        from scipy.stats import poisson
        x = np.array([counts.get(s, 0) for s in self.species], dtype=np.int64)
        if not self.is_conditioned:
            return float(np.prod([poisson.pmf(xi, li)
                                  for xi, li in zip(x, self.means)]))
        if self._states is None:
            # Could not enumerate the class → independent-Poisson approximation.
            return float(np.prod([poisson.pmf(xi, li)
                                  for xi, li in zip(x, self.means)]))
        from scipy.special import logsumexp
        if self._logZ is None:
            self._logZ = float(logsumexp([self._log_weight(s) for s in self._states]))
        return float(np.exp(self._log_weight(x) - self._logZ))

    def expected_counts(self) -> dict[str, float]:
        """Stationary mean molecule counts (conditional means when conditioned)."""
        if not self.is_conditioned or self._states is None:
            return {s: float(m) for s, m in zip(self.species, self.means)}
        probs = np.array([self.probability(dict(zip(self.species, st)))
                          for st in self._states])
        mean = (self._states * probs[:, None]).sum(axis=0)
        return {s: float(m) for s, m in zip(self.species, mean)}

    def marginal(self, species: str) -> tuple[np.ndarray, np.ndarray]:
        """Return ``(values, probabilities)`` for the marginal of one species."""
        i = self.species.index(species)
        if not self.is_conditioned or self._states is None:
            from scipy.stats import poisson
            hi = int(poisson.ppf(0.9999, self.means[i])) + 1
            vals = np.arange(0, hi + 1)
            return vals, poisson.pmf(vals, self.means[i])
        probs = np.array([self.probability(dict(zip(self.species, st)))
                          for st in self._states])
        col = self._states[:, i]
        vals = np.arange(col.min(), col.max() + 1)
        out = np.zeros(len(vals))
        for v, p in zip(col, probs):
            out[v - col.min()] += p
        return vals, out

    def sample(self, size: int = 1, seed: int | None = None) -> dict[str, np.ndarray]:
        """Draw samples from the stationary distribution."""
        rng = np.random.default_rng(seed)
        if not self.is_conditioned or self._states is None:
            draws = np.array([rng.poisson(m, size=size) for m in self.means]).T
        else:
            probs = np.array([self.probability(dict(zip(self.species, st)))
                              for st in self._states])
            probs = probs / probs.sum()
            idx = rng.choice(len(self._states), size=size, p=probs)
            draws = self._states[idx]
        return {s: draws[:, i] for i, s in enumerate(self.species)}

expected_counts()

Stationary mean molecule counts (conditional means when conditioned).

Source code in mantis/stochastic_stationary.py

def expected_counts(self) -> dict[str, float]:
    """Stationary mean molecule counts (conditional means when conditioned)."""
    if not self.is_conditioned or self._states is None:
        return {s: float(m) for s, m in zip(self.species, self.means)}
    probs = np.array([self.probability(dict(zip(self.species, st)))
                      for st in self._states])
    mean = (self._states * probs[:, None]).sum(axis=0)
    return {s: float(m) for s, m in zip(self.species, mean)}

marginal(species)

Return (values, probabilities) for the marginal of one species.

Source code in mantis/stochastic_stationary.py

def marginal(self, species: str) -> tuple[np.ndarray, np.ndarray]:
    """Return ``(values, probabilities)`` for the marginal of one species."""
    i = self.species.index(species)
    if not self.is_conditioned or self._states is None:
        from scipy.stats import poisson
        hi = int(poisson.ppf(0.9999, self.means[i])) + 1
        vals = np.arange(0, hi + 1)
        return vals, poisson.pmf(vals, self.means[i])
    probs = np.array([self.probability(dict(zip(self.species, st)))
                      for st in self._states])
    col = self._states[:, i]
    vals = np.arange(col.min(), col.max() + 1)
    out = np.zeros(len(vals))
    for v, p in zip(col, probs):
        out[v - col.min()] += p
    return vals, out

poisson_means()

Poisson means λ_i keyed by species (molecule counts).

Source code in mantis/stochastic_stationary.py

def poisson_means(self) -> dict[str, float]:
    """Poisson means ``λ_i`` keyed by species (molecule counts)."""
    return {s: float(m) for s, m in zip(self.species, self.means)}

probability(counts)

Normalised stationary probability of an integer state counts.

Source code in mantis/stochastic_stationary.py

def probability(self, counts: dict[str, float]) -> float:
    """Normalised stationary probability of an integer state ``counts``."""
    from scipy.stats import poisson
    x = np.array([counts.get(s, 0) for s in self.species], dtype=np.int64)
    if not self.is_conditioned:
        return float(np.prod([poisson.pmf(xi, li)
                              for xi, li in zip(x, self.means)]))
    if self._states is None:
        # Could not enumerate the class → independent-Poisson approximation.
        return float(np.prod([poisson.pmf(xi, li)
                              for xi, li in zip(x, self.means)]))
    from scipy.special import logsumexp
    if self._logZ is None:
        self._logZ = float(logsumexp([self._log_weight(s) for s in self._states]))
    return float(np.exp(self._log_weight(x) - self._logZ))

sample(size=1, seed=None)

Draw samples from the stationary distribution.

Source code in mantis/stochastic_stationary.py

def sample(self, size: int = 1, seed: int | None = None) -> dict[str, np.ndarray]:
    """Draw samples from the stationary distribution."""
    rng = np.random.default_rng(seed)
    if not self.is_conditioned or self._states is None:
        draws = np.array([rng.poisson(m, size=size) for m in self.means]).T
    else:
        probs = np.array([self.probability(dict(zip(self.species, st)))
                          for st in self._states])
        probs = probs / probs.sum()
        idx = rng.choice(len(self._states), size=size, p=probs)
        draws = self._states[idx]
    return {s: draws[:, i] for i, s in enumerate(self.species)}

stationary_distribution(reactions, species, rate_values, deficiency, is_weakly_reversible, initial_conditions, volume_L, initial_as='concentration', chemostatted_values=None, max_states=200000)

Build the exact ACK product-form stationary distribution.

Raises ValueError if the network is not complex-balanced (the closed-form result does not apply).

Source code in mantis/stochastic_stationary.py

def stationary_distribution(
    reactions: list[Reaction],
    species: list[str],
    rate_values: dict[str, float],
    deficiency: int,
    is_weakly_reversible: bool,
    initial_conditions: dict[str, float],
    volume_L: float,
    initial_as: str = "concentration",
    chemostatted_values: dict[str, float] | None = None,
    max_states: int = 200_000,
) -> StationaryDistribution:
    """Build the exact ACK product-form stationary distribution.

    Raises ``ValueError`` if the network is not complex-balanced (the closed-form
    result does not apply).
    """
    from .stability import is_complex_balanced, complex_balanced_equilibrium
    from .analysis import _conservation_law_vectors

    chem = chemostatted_values or {}
    NA_V = AVOGADRO * volume_L

    # Initial molecule counts and the equivalent concentrations.
    n0 = np.zeros(len(species))
    for i, s in enumerate(species):
        val = initial_conditions.get(s, 0.0)
        n0[i] = round(val) if initial_as == "count" else round(val * NA_V)
    ic_conc = {s: n0[i] / NA_V for i, s in enumerate(species)}

    if not is_complex_balanced(reactions, species, rate_values, deficiency,
                               is_weakly_reversible, ic_conc, chem or None):
        raise ValueError(
            "Network is not complex-balanced at these rate constants, so the "
            "Anderson–Craciun–Kurtz product-form stationary distribution does not "
            "apply. (It holds for every weakly reversible deficiency-zero network, "
            "and for higher-deficiency networks only at complex-balanced rates.)"
        )

    c_star = complex_balanced_equilibrium(
        reactions, species, rate_values, deficiency, is_weakly_reversible,
        ic_conc, chem or None,
    )
    means = np.array([c_star[s] * NA_V for s in species])

    cons_vectors = _conservation_law_vectors(reactions, species)
    is_conditioned = len(cons_vectors) > 0
    states = None
    is_exact = True
    if is_conditioned:
        cons_totals = [float(v @ n0) for v in cons_vectors]
        states, exact = _enumerate_class(cons_vectors, cons_totals, len(species), max_states)
        is_exact = exact

    return StationaryDistribution(
        species=list(species),
        means=means,
        volume_L=volume_L,
        is_conditioned=is_conditioned,
        is_exact=is_exact,
        _states=states,
    )

Finite State Projection

mantis.fsp

Finite State Projection — Munsky & Khammash (J. Chem. Phys. 2006).

The chemical master equation (CME) is an (often infinite) linear ODE for the probability of every molecular-count state. FSP truncates the state space to a finite, reachable subset, assembles the corresponding generator matrix, and solves the truncated CME directly — giving accurate distributions exactly where the SSA is noisiest (the low-copy-number biosensor regime), together with a rigorous bound on the truncation error (the probability that has leaked out of the projection).

This module enumerates the reachable states (BFS under the reaction stoichiometry, respecting conservation laws, capped at max_states), builds the sparse generator Q, and provides:

• transient distributions p(t) = expm(Q·t)·p₀ with the leaked-mass error bound;
• stationary distributions (null space of Q) for closed, fully enumerated networks — which, for complex-balanced networks, reproduces the exact Anderson–Craciun–Kurtz product form (a strong cross-check).

FSPResult dataclass

A finite-state-projection distribution over molecular-count states.

Attributes

species : list[str] states : np.ndarray (M, n) integer array of enumerated states. probabilities : np.ndarray Length-M probability vector aligned with states. truncation_error : float Probability that left the projection (0 for a closed, fully enumerated network; an upper bound on the distribution error otherwise). time : float | None Evaluation time for a transient solve, or None for the stationary one. truncated : bool True if state enumeration hit max_states before closing.

Source code in mantis/fsp.py

@dataclass
class FSPResult:
    """A finite-state-projection distribution over molecular-count states.

    Attributes
    ----------
    species : list[str]
    states : np.ndarray
        ``(M, n)`` integer array of enumerated states.
    probabilities : np.ndarray
        Length-``M`` probability vector aligned with ``states``.
    truncation_error : float
        Probability that left the projection (0 for a closed, fully enumerated
        network; an upper bound on the distribution error otherwise).
    time : float | None
        Evaluation time for a transient solve, or ``None`` for the stationary one.
    truncated : bool
        True if state enumeration hit ``max_states`` before closing.
    """
    species: list[str]
    states: np.ndarray
    probabilities: np.ndarray
    truncation_error: float
    time: float | None = None
    truncated: bool = False
    _index: dict = field(default=None, repr=False)

    def probability(self, counts: dict[str, float]) -> float:
        """Probability of an exact integer state (0 if outside the projection)."""
        key = tuple(int(counts.get(s, 0)) for s in self.species)
        if self._index is None:
            self._index = {tuple(int(v) for v in row): i
                           for i, row in enumerate(self.states)}
        i = self._index.get(key)
        return float(self.probabilities[i]) if i is not None else 0.0

    def marginal(self, species: str) -> tuple[np.ndarray, np.ndarray]:
        """Return ``(values, probabilities)`` for one species' marginal."""
        i = self.species.index(species)
        col = self.states[:, i]
        vals = np.arange(col.min(), col.max() + 1)
        out = np.zeros(len(vals))
        for v, p in zip(col, self.probabilities):
            out[v - col.min()] += p
        return vals, out

    def expected_counts(self) -> dict[str, float]:
        """Mean molecule counts under the (renormalised) projected distribution."""
        w = self.probabilities / self.probabilities.sum()
        mean = (self.states * w[:, None]).sum(axis=0)
        return {s: float(m) for s, m in zip(self.species, mean)}

expected_counts()

Mean molecule counts under the (renormalised) projected distribution.

Source code in mantis/fsp.py

def expected_counts(self) -> dict[str, float]:
    """Mean molecule counts under the (renormalised) projected distribution."""
    w = self.probabilities / self.probabilities.sum()
    mean = (self.states * w[:, None]).sum(axis=0)
    return {s: float(m) for s, m in zip(self.species, mean)}

marginal(species)

Return (values, probabilities) for one species' marginal.

Source code in mantis/fsp.py

def marginal(self, species: str) -> tuple[np.ndarray, np.ndarray]:
    """Return ``(values, probabilities)`` for one species' marginal."""
    i = self.species.index(species)
    col = self.states[:, i]
    vals = np.arange(col.min(), col.max() + 1)
    out = np.zeros(len(vals))
    for v, p in zip(col, self.probabilities):
        out[v - col.min()] += p
    return vals, out

probability(counts)

Probability of an exact integer state (0 if outside the projection).

Source code in mantis/fsp.py

def probability(self, counts: dict[str, float]) -> float:
    """Probability of an exact integer state (0 if outside the projection)."""
    key = tuple(int(counts.get(s, 0)) for s in self.species)
    if self._index is None:
        self._index = {tuple(int(v) for v in row): i
                       for i, row in enumerate(self.states)}
    i = self._index.get(key)
    return float(self.probabilities[i]) if i is not None else 0.0

fsp_solve(reactions, species, rate_values, initial_conditions, volume_L, t=None, initial_as='concentration', chemostatted_values=None, max_states=100000)

Solve the CME by finite state projection.

With t given, returns the transient distribution p(t); with t=None, returns the stationary distribution (null space of Q) — valid when the enumerated state space is closed.

Source code in mantis/fsp.py

def fsp_solve(
    reactions: list[Reaction],
    species: list[str],
    rate_values: dict[str, float],
    initial_conditions: dict[str, float],
    volume_L: float,
    t: float | None = None,
    initial_as: str = "concentration",
    chemostatted_values: dict[str, float] | None = None,
    max_states: int = 100_000,
) -> FSPResult:
    """Solve the CME by finite state projection.

    With ``t`` given, returns the transient distribution ``p(t)``; with ``t=None``,
    returns the stationary distribution (null space of ``Q``) — valid when the
    enumerated state space is closed.
    """
    na_v = AVOGADRO * volume_L
    n0 = np.zeros(len(species), dtype=np.int64)
    for i, s in enumerate(species):
        val = initial_conditions.get(s, 0.0)
        n0[i] = int(round(val)) if initial_as == "count" else int(round(val * na_v))

    c_eff, react_sp, react_co, change = _kernel_arrays(
        reactions, species, rate_values, chemostatted_values, volume_L
    )
    propensities = _propensity_fn(c_eff, react_sp, react_co)
    states, index, truncated = _enumerate_reachable(n0, change, propensities, max_states)
    Q, leak = _build_generator(states, index, change, propensities)
    M = len(states)

    if t is not None:
        from scipy.sparse.linalg import expm_multiply
        p0 = np.zeros(M)
        p0[index[tuple(int(v) for v in n0)]] = 1.0
        p = expm_multiply(Q * float(t), p0)
        p = np.maximum(p, 0.0)
        trunc_err = float(max(0.0, 1.0 - p.sum()))
        return FSPResult(
            species=list(species), states=states, probabilities=p,
            truncation_error=trunc_err, time=float(t), truncated=truncated,
            _index=index,
        )

    # Stationary distribution: solve Qπ = 0 with Σπ = 1 (replace one row with the
    # normalisation constraint).
    A = Q.toarray()
    A[-1, :] = 1.0
    b = np.zeros(M)
    b[-1] = 1.0
    pi = np.linalg.solve(A, b)
    pi = np.maximum(pi, 0.0)
    pi = pi / pi.sum()
    trunc_err = float(leak.max()) if M else 0.0
    return FSPResult(
        species=list(species), states=states, probabilities=pi,
        truncation_error=trunc_err, time=None, truncated=truncated,
        _index=index,
    )

SBML interoperability

mantis.sbml

SBML (Systems Biology Markup Language) import/export for CRNetwork.

This is an optional interoperability layer. It requires python-libsbml; if that package is not installed, importing the helpers raises a clear :class:ImportError with an install hint, but the rest of mantis keeps working without it.

The mapping between mantis and SBML is mass-action-centric:

• Each directed :class:~mantis.parsing.Reaction becomes one irreversible SBML reaction. A reversible A <-> B round-trips as two irreversible reactions (matching mantis' internal representation).
• Each reaction carries a mass-action kinetic law k_i * Π reactant^coeff with the rate constant stored as a global parameter k_i.
• Chemostatted species are exported as SBML boundary species (boundaryCondition=True) holding their fixed concentration; on import, boundary / constant species are recovered as chemostatted.

network_from_sbml(source)

Construct a :class:CRNetwork from an SBML file path or XML string.

source is treated as a file path if it points to an existing file, otherwise as a raw SBML/XML document string.

Source code in mantis/sbml.py

def network_from_sbml(source: str) -> "CRNetwork":
    """Construct a :class:`CRNetwork` from an SBML file path or XML string.

    ``source`` is treated as a file path if it points to an existing file,
    otherwise as a raw SBML/XML document string.
    """
    import os

    from .network import CRNetwork
    from .parsing import parse_reactions, normalize_rate_key

    libsbml = _require_libsbml()

    if os.path.exists(source):
        doc = libsbml.readSBMLFromFile(source)
    else:
        doc = libsbml.readSBMLFromString(source)
    if doc is None or doc.getModel() is None:
        raise ValueError("Could not parse SBML: no model found in document.")
    if doc.getNumErrors(libsbml.LIBSBML_SEV_ERROR) > 0:
        msg = doc.getErrorLog().toString()
        raise ValueError(f"SBML document has fatal errors:\n{msg}")

    model = doc.getModel()
    species_ids = {model.getSpecies(i).getId() for i in range(model.getNumSpecies())}

    # Chemostatted species: SBML boundary or constant species, with their fixed value.
    chemostatted: dict[str, float] = {}
    for i in range(model.getNumSpecies()):
        sp = model.getSpecies(i)
        if sp.getBoundaryCondition() or sp.getConstant():
            if sp.isSetInitialConcentration():
                val = float(sp.getInitialConcentration())
            elif sp.isSetInitialAmount():
                val = float(sp.getInitialAmount())
            else:
                val = 0.0
            chemostatted[sp.getId()] = val

    reaction_strings: list[str] = []
    rates: dict[str, float] = {}
    for i in range(model.getNumReactions()):
        rxn = model.getReaction(i)
        reactants = _stoich_terms(
            [rxn.getReactant(j) for j in range(rxn.getNumReactants())]
        )
        products = _stoich_terms(
            [rxn.getProduct(j) for j in range(rxn.getNumProducts())]
        )
        lhs = _complex_str(reactants)
        rhs = _complex_str(products)
        arrow = "<->" if rxn.getReversible() else "->"
        rxn_str = f"{lhs} {arrow} {rhs}"
        reaction_strings.append(rxn_str)

        fwd, rev = _rates_for_reaction(model, rxn, species_ids, libsbml)
        if fwd is not None:
            rates[normalize_rate_key(f"{lhs} -> {rhs}")] = fwd
        if rev is not None:
            rates[normalize_rate_key(f"{rhs} -> {lhs}")] = rev

    reactions = parse_reactions(reaction_strings)
    return CRNetwork(reactions, rates=rates, chemostatted=chemostatted or None)

network_to_sbml(net, *, level=3, version=2)

Serialize a :class:CRNetwork to an SBML document string (Level 3 v2).

Source code in mantis/sbml.py

def network_to_sbml(net: "CRNetwork", *, level: int = 3, version: int = 2) -> str:
    """Serialize a :class:`CRNetwork` to an SBML document string (Level 3 v2)."""
    libsbml = _require_libsbml()

    doc = libsbml.SBMLDocument(level, version)
    model = doc.createModel()
    model.setId("mantis_model")

    # A single well-mixed compartment of unit size.
    comp = model.createCompartment()
    comp.setId("c")
    comp.setConstant(True)
    comp.setSize(1.0)
    comp.setSpatialDimensions(3)

    chem = net.chemostatted  # name -> fixed concentration

    # All species names appearing anywhere in the reactions (dynamic + chemostatted).
    names: list[str] = []
    seen: set[str] = set()
    for rxn in net._reactions:
        for name, _ in list(rxn.reactants) + list(rxn.products):
            if name not in seen:
                seen.add(name)
                names.append(name)
    names.sort()

    for name in names:
        sp = model.createSpecies()
        sp.setId(name)
        sp.setCompartment("c")
        sp.setHasOnlySubstanceUnits(False)
        if name in chem:
            sp.setInitialConcentration(float(chem[name]))
            sp.setBoundaryCondition(True)
            sp.setConstant(True)
        else:
            sp.setInitialConcentration(0.0)
            sp.setBoundaryCondition(False)
            sp.setConstant(False)

    # One global rate parameter per reaction, plus one irreversible reaction.
    for i, rxn in enumerate(net._reactions, start=1):
        kid = f"k_{i}"
        param = model.createParameter()
        param.setId(kid)
        param.setConstant(True)
        # Only assign a value for rates that were actually supplied; leaving it
        # unset lets ``from_sbml`` distinguish "unknown" from a literal 1.0.
        if rxn.rate_key in net._rates:
            param.setValue(float(net._rates[rxn.rate_key]))

        r = model.createReaction()
        r.setId(f"r_{i}")
        r.setReversible(False)
        # The reaction's human-readable rate key is preserved as the SBML name.
        r.setName(rxn.rate_key)

        for name, coeff in sorted(rxn.reactants, key=lambda x: x[0]):
            sr = r.createReactant()
            sr.setSpecies(name)
            sr.setStoichiometry(float(coeff))
            sr.setConstant(True)
        for name, coeff in sorted(rxn.products, key=lambda x: x[0]):
            sp_ref = r.createProduct()
            sp_ref.setSpecies(name)
            sp_ref.setStoichiometry(float(coeff))
            sp_ref.setConstant(True)

        # Mass-action kinetic law: k_i * Π species^coeff (chemostatted included).
        factors = [kid]
        for name, coeff in sorted(rxn.reactants, key=lambda x: x[0]):
            factors.extend([name] * int(coeff))
        formula = " * ".join(factors)
        kl = r.createKineticLaw()
        math_ast = libsbml.parseL3Formula(formula)
        if math_ast is None:  # pragma: no cover - formula is machine-generated
            raise RuntimeError(f"Failed to parse kinetic-law formula: {formula!r}")
        kl.setMath(math_ast)

    writer = libsbml.SBMLWriter()
    return writer.writeSBMLToString(doc)

write_sbml(net, path)

Write a :class:CRNetwork to an SBML file at path.

Source code in mantis/sbml.py

def write_sbml(net: "CRNetwork", path: str) -> None:
    """Write a :class:`CRNetwork` to an SBML file at ``path``."""
    xml = network_to_sbml(net)
    with open(path, "w", encoding="utf-8") as fh:
        fh.write(xml)