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 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 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 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 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)

chemostatted property

Chemostatted species and their fixed concentrations.

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

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)

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)

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"]

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,
    )

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,
    )

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

    rng = np.random.default_rng(seed)
    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
    n = 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":
            n[sp_idx[s]] = int(round(val))
        else:
            n[sp_idx[s]] = int(round(val * NA_V))

    # Per-reaction reactant info and stoichiometric change vector
    rxn_info = []  # list of (k_deterministic, [(idx, coeff), ...], change_vec)
    for rxn in reactions:
        k_det = float(eff_rates.get(rxn.rate_key, 0.0))
        reactants = [(sp_idx[name], coeff) for name, coeff in rxn.reactants
                     if name not in chem]
        change = np.zeros(n_sp, dtype=np.int64)
        for name, coeff in rxn.reactants:
            if name in sp_idx and name not in chem:
                change[sp_idx[name]] -= coeff
        for name, coeff in rxn.products:
            if name in sp_idx and name not in chem:
                change[sp_idx[name]] += coeff
        rxn_info.append((k_det, reactants, change))

    def propensity(rxn_idx: int, state: np.ndarray) -> float:
        """Stochastic propensity for reaction rxn_idx given current counts."""
        k_det, reactants, _ = rxn_info[rxn_idx]
        if k_det <= 0.0:
            return 0.0
        order = sum(c for _, c in reactants)
        # h(x) = product of combinatorial choose terms
        h = 1.0
        for idx, coeff in reactants:
            ni = int(state[idx])
            if ni < coeff:
                return 0.0
            for j in range(coeff):
                h *= (ni - j)
            if coeff > 1:
                h /= math.factorial(coeff)
        # c = k_det / (N_A V)^(order - 1)
        if order > 1:
            c = k_det / (NA_V ** (order - 1))
        else:
            c = k_det
        return c * h

    times = [float(t_span[0])]
    history = [n.copy()]
    t = float(t_span[0])
    tf = float(t_span[1])
    n_events = 0
    success = True

    # Pre-allocate propensity array
    a = np.zeros(len(reactions))

    while t < tf:
        if n_events >= max_events:
            success = False
            break
        for j in range(len(reactions)):
            a[j] = propensity(j, n)
        a0 = float(a.sum())
        if a0 <= 0.0:
            # No more reactions can fire; system frozen.
            break
        # Time to next event
        r1, r2 = rng.random(), rng.random()
        tau = -math.log(r1) / a0
        t_new = t + tau
        if t_new > tf:
            break
        # Pick reaction
        threshold = r2 * a0
        cum = 0.0
        chosen = len(reactions) - 1
        for j in range(len(reactions)):
            cum += a[j]
            if cum >= threshold:
                chosen = j
                break
        # Apply
        n = n + rxn_info[chosen][2]
        t = t_new
        times.append(t)
        history.append(n.copy())
        n_events += 1

    # Pad to tf
    if times[-1] < tf:
        times.append(tf)
        history.append(n.copy())

    times_arr = np.asarray(times)
    history_arr = np.stack(history, axis=0)   # shape (T, n_sp)
    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=n_events,
        success=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

    rng = np.random.default_rng(seed)
    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

    n = 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":
            n[sp_idx[s]] = int(round(val))
        else:
            n[sp_idx[s]] = int(round(val * NA_V))

    rxn_info = []
    change_matrix = np.zeros((len(reactions), n_sp), dtype=np.int64)
    for ridx, rxn in enumerate(reactions):
        k_det = float(eff_rates.get(rxn.rate_key, 0.0))
        reactants = [(sp_idx[name], coeff) for name, coeff in rxn.reactants
                     if name not in chem]
        for name, coeff in rxn.reactants:
            if name in sp_idx and name not in chem:
                change_matrix[ridx, sp_idx[name]] -= coeff
        for name, coeff in rxn.products:
            if name in sp_idx and name not in chem:
                change_matrix[ridx, sp_idx[name]] += coeff
        rxn_info.append((k_det, reactants))

    def propensity(rxn_idx: int, state: np.ndarray) -> float:
        k_det, reactants = rxn_info[rxn_idx]
        if k_det <= 0.0:
            return 0.0
        order = sum(c for _, c in reactants)
        h = 1.0
        for idx, coeff in reactants:
            ni = int(state[idx])
            if ni < coeff:
                return 0.0
            for j in range(coeff):
                h *= (ni - j)
            if coeff > 1:
                h /= math.factorial(coeff)
        if order > 1:
            c = k_det / (NA_V ** (order - 1))
        else:
            c = k_det
        return c * h

    def all_propensities(state: np.ndarray) -> np.ndarray:
        out = np.empty(len(reactions))
        for j in range(len(reactions)):
            out[j] = propensity(j, state)
        return out

    def adaptive_tau(state: np.ndarray, a: np.ndarray) -> float:
        """Cao et al. 2006: bound the relative change of each propensity."""
        a0 = float(a.sum())
        if a0 <= 0:
            return float("inf")
        # Per-species drift and variance under one leap
        # mu_i = Σ_j ν_ji · a_j ;  σ²_i = Σ_j ν_ji² · a_j
        nu = change_matrix  # (n_rxn, n_sp)
        mu = nu.T @ a       # (n_sp,)
        sigma2 = (nu.T ** 2) @ a
        # Highest-order reaction touching each species (Cao approximation g_i)
        # Simplification: g_i = 2 for all species (works for ≤ 2nd-order rxns)
        g = np.full(n_sp, 2.0)
        x = np.maximum(state.astype(float), 1.0)
        bound1 = np.where(np.abs(mu) > 0, np.maximum(epsilon * x / g, 1.0) / np.abs(mu), np.inf)
        bound2 = np.where(sigma2 > 0, (np.maximum(epsilon * x / g, 1.0) ** 2) / sigma2, np.inf)
        return float(min(bound1.min(), bound2.min()))

    t = float(t_span[0])
    tf = float(t_span[1])
    record_times = np.linspace(t, tf, max(2, n_record))
    record_states = np.zeros((len(record_times), n_sp), dtype=np.int64)
    record_states[0] = n
    next_record = 1
    success = True

    while t < tf and next_record < len(record_times):
        a = all_propensities(n)
        a0 = float(a.sum())
        if a0 <= 0.0:
            break

        if tau is not None:
            dt = float(tau)
        else:
            dt = adaptive_tau(n, a)

        # Cap the step to the next recording time so we don't skip past it.
        dt = min(dt, tf - t, record_times[next_record] - t)
        if dt <= 0:
            dt = tf - t

        # Sample firing counts
        firings = rng.poisson(a * dt)

        # Provisional new state — reject if it goes negative
        delta = firings @ change_matrix
        n_new = n + delta
        if (n_new < 0).any():
            # Leap rejection: fall back to one exact SSA step
            r1, r2 = rng.random(), rng.random()
            ssa_dt = -math.log(max(r1, 1e-300)) / a0
            if t + ssa_dt > tf:
                t = tf
                break
            cum = 0.0
            threshold = r2 * a0
            chosen = len(reactions) - 1
            for j in range(len(reactions)):
                cum += a[j]
                if cum >= threshold:
                    chosen = j
                    break
            n = n + change_matrix[chosen]
            t = t + ssa_dt
        else:
            n = n_new
            t = t + dt

        # Snapshot at any recording times we've crossed
        while next_record < len(record_times) and record_times[next_record] <= t:
            record_states[next_record] = n
            next_record += 1

    # Fill any remaining slots with the last state
    for i in range(next_record, len(record_times)):
        record_states[i] = n

    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,
    )