Skip to content

iocrn

RL4CRN.iocrns.iocrn

Input-Output Chemical Reaction Network (IOCRN) representation.

This module defines IOCRN, a representation of an input-output chemical reaction network with utilities for:

  • building and compiling a CRN from a list of reactions,
  • mapping between symbolic labels (species/inputs) and numeric indices,
  • simulating transient responses (ODE solvers and SSA),
  • plotting simulation results,
  • exporting the network to a DSL format used by external simulators,
  • (optionally) solving dynamics using CVODE via sksundae.

This modules provides utilities for deterministic simulations of CRNs, while it relies on an external package for stochastic simulation (SSA).

Caching

Simulation results are cached in last_task_info to avoid recomputation if the same task is requested repeatedly. Adding reactions or calling reset clears cached task info.

Solver support
  • 'LSODA' uses scipy.integrate.solve_ivp.
  • 'CVODE' uses sksundae.cvode.CVODE.
SSA support

Stochastic simulation is supported via an external package StochasticSimulationsNew plus project utilities in :mod:RL4CRN.utils.stochastic. If the external package is missing, SSA methods raise an ImportError.

IOCRN

Input-Output Chemical Reaction Network.

An IOCRN consists of
  • a list of reactions,
  • a set of output species labels,
  • an implied set of species and inputs gathered from reactions,
  • dynamics defined by stoichiometry and propensities.

The network must be compiled after construction or modification to build internal indices and matrices used for simulation.

PARAMETER DESCRIPTION
reactions

List of Reaction objects. Each reaction is expected to provide methods such as:

  • get_involved_species()
  • get_involved_inputs()
  • get_stoichiometry_dict()
  • propensity(x, u)
  • set_crn_context(iocrn)

and attributes such as:

  • ID
  • params
  • num_parameters
  • num_unknown_params
output_labels

List of species labels treated as outputs.
solver

ODE solver identifier. Supported values:

  • 'CVODE' (requires sksundae)
  • 'LSODA' (SciPy solve_ivp)

DEFAULT: 'CVODE'

atol

Absolute tolerance for ODE integration.

DEFAULT: 1e-06

rtol

Relative tolerance for ODE integration.

DEFAULT: 0.001

ATTRIBUTE DESCRIPTION
reactions

List of reaction objects.
output_labels

List of output species labels.
num_outputs

Number of outputs.
num_unknown_params

Total number of unknown parameters across reactions.
last_task_info

Dictionary caching results/metadata of the last evaluated task (e.g., transient response).
solver

Selected solver backend name.
atol

Absolute tolerance.
rtol

Relative tolerance.

__init__(reactions, output_labels, solver='CVODE', atol=1e-06, rtol=0.001)

Initialize the IOCRN.

Important

After initialization, call compile to create internal representations (species/input indices, stoichiometry matrix, etc.).

PARAMETER DESCRIPTION
reactions

List of reaction objects.
output_labels

List of labels for output species.
solver

ODE solver backend ('CVODE' or 'LSODA').

DEFAULT: 'CVODE'

atol

Absolute tolerance passed to the solver.

DEFAULT: 1e-06

rtol

Relative tolerance passed to the solver.

DEFAULT: 0.001

clone()

Return a deep copy of the IOCRN.

reset()

Clear cached task information.

This does not modify reactions or parameters; it only clears the cached last_task_info.

add_reaction(reaction)

Add a reaction and recompile.

Adding a reaction clears cached task info and updates internal representations by calling compile.

PARAMETER DESCRIPTION
reaction

Reaction object to append to the network.

retrieve_species()

Collect species labels from reactions and create index mappings.

Side Effects

Sets:

  • species_labels: sorted list of all species labels
  • species_idx_dict: mapping label -> index

retrieve_input()

Collect input labels from reactions and create index mappings.

Side Effects

Sets:

  • input_labels: sorted list of all input labels
  • num_inputs: number of inputs
  • input_idx_dict: mapping label -> index

species_label_to_idx(labels)

Map species labels to indices.

PARAMETER DESCRIPTION
labels

Species label (str) or list of species labels.
RETURNS DESCRIPTION

If labels is a string, returns an integer index. If labels is a list, returns a list of integer indices.

input_label_to_idx(labels)

Map input labels to indices.

PARAMETER DESCRIPTION
labels

Input label (str) or list of input labels. Elements may be None, which are preserved.
RETURNS DESCRIPTION

If labels is a string, returns an integer index. If labels is a list, returns a list of indices (with None preserved).

set_library_context(reaction_library)

Set a reaction-library context for all reactions.

This is required for some topology comparisons and signatures.

PARAMETER DESCRIPTION
reaction_library

Library object passed to reaction.set_library_context(...).
Side Effects

Sets self.reaction_library.

get_bool_signature()

Return a boolean signature of present reactions relative to a library.

RETURNS DESCRIPTION

Boolean numpy array of length len(self.reaction_library) where True indicates the corresponding library reaction ID is present.
RAISES DESCRIPTION
AttributeError

If reaction_library has not been set.

is_topologically_equal(other_iocrn)

Check topology equality against another IOCRN.

Assumes both IOCRNs share the same reaction library context.

PARAMETER DESCRIPTION
other_iocrn

Another IOCRN.
RETURNS DESCRIPTION

True if their boolean signatures match, else False.

get_stoichiometry_matrix()

Construct the stoichiometry matrix \(S\).

RETURNS DESCRIPTION

Numpy array S of shape (num_species, num_reactions) where S[i, j] is the net stoichiometric coefficient of species i in reaction j.

gather_reaction_IDs()

Return reaction IDs for all reactions in the network.

gather_reaction_params()

Return a list of parameter vectors for all reactions.

compile()

Compile internal IOCRN representations.

This method
  • retrieves and indexes species and inputs,
  • computes output_idx for output species,
  • builds the stoichiometry matrix S,
  • calls reaction.set_crn_context(self) on each reaction.
Call this after
  • initialization,
  • adding/removing reactions,
  • changing reaction definitions that affect involved species/inputs.
Important

This method should be called after adding all reactions and before simulating the IOCRN.

__str__()

Return a readable string representation of the IOCRN.

Includes
  • inputs
  • species
  • output species
  • reactions (sorted by reaction ID if available)

to_reaction_file()

Export the IOCRN to a custom DSL reaction-file format.

The exported text contains
  1. An input ...; section (if inputs are present),
  2. A species ... = 0; section (excluding emptyset/∅),
  3. A reactions section where each reaction contributes one or more lines via reaction.to_reaction_format().
RETURNS DESCRIPTION

The DSL text as a single string.

propensity_function(x, u)

Compute propensities \(a(x,u)\) for all reactions.

PARAMETER DESCRIPTION
x

Species state vector (array-like of shape (n,)).
u

Input vector (array-like of shape (p,)).
RETURNS DESCRIPTION

Numpy array of propensities of shape (num_reactions,).

rate_function(t, x, u)

Compute time-derivative \(\dot{x} = S a(x,u)\).

PARAMETER DESCRIPTION
t

Time (unused for autonomous dynamics, but accepted by solver APIs).
x

Species state vector of shape (n,).
u

Input vector of shape (p,).
RETURNS DESCRIPTION

Numpy array of shape (n,) giving \(\dot{x}\).

dose_response(u_dose, u_list, initial_guess)

Computes the dose response of the IOCRN given a list of input doses, a list of input scenarios, and an initial guess for the concentrations. If the CRN dose response has been simulated and stored before, it returns the stored results instead of recomputing them. The results are stored in the last_task_info dictionary for future reference. Arguments: - u_dose: numpy array of shape (num_doses,) representing the input doses to the IOCRN. - u_list: A list of numpy arrays, each of shape (p,) representing the constant inputs to the IOCRN for each scenario. The first element of each input array corresponds to the dose, and the rest correspond to other inputs. - initial_guess: numpy array of shape (n,) representing the initial guess for the concentrations of the species. Returns a tupple containing: - x_list: A list of numpy arrays of shape (n, num_doses) representing the state for each input scenario. - y_list: A list of numpy arrays of shape (q, num_doses) representing the output for each input scenario.

transient_response_SSA(u_list, x0_list, time_horizon, n_trajectories=100, max_threads=10000, max_value=1000000.0)

Compute stochastic transient responses via SSA (mean ± std).

This method runs a stochastic simulation algorithm (SSA) backend for the current IOCRN by first exporting the network to the project DSL (to_reaction_file) and then invoking the external SSA engine.

For each pair (u, x0) in the Cartesian product u_list × x0_list, the method performs n_trajectories independent SSA runs and returns the empirical mean and standard deviation trajectories for all species and outputs. Results are interpolated onto the requested time_horizon.

The results are cached in last_task_info with type == 'transient response SSA'.

PARAMETER DESCRIPTION
u_list

List of input vectors. Each element is array-like of shape (p,), where p = num_inputs. Inputs are treated as constant over time for each scenario.
x0_list

List of initial condition vectors. Each element is array-like of shape (n,), where n = num_species.
time_horizon

1D array-like of time points at which to return statistics. The SSA backend is sampled at a step t_step inferred from time_horizon and then interpolated to match it.
n_trajectories

Number of SSA trajectories per (u, x0) scenario used to estimate mean and standard deviation.

DEFAULT: 100

max_threads

Maximum threads/parallelism hint for the SSA backend. (Currently not directly used in this method; kept for API compatibility / future backend control.)

DEFAULT: 10000

max_value

Divergence threshold passed to the SSA measurement helper. Trajectories exceeding this value may be flagged as diverged.

DEFAULT: 1000000.0

RETURNS DESCRIPTION

Tuple (time_horizon, x_mean_list, y_mean_list, x_std_list, y_std_list, last_task_info):

  • time_horizon: The same array of time points passed in.
  • x_mean_list: List of mean state trajectories, one per scenario. Each entry has shape (n, T).
  • y_mean_list: List of mean output trajectories, one per scenario. Each entry has shape (q, T).
  • x_std_list: List of std-dev state trajectories, one per scenario. Each entry has shape (n, T).
  • y_std_list: List of std-dev output trajectories, one per scenario. Each entry has shape (q, T).
  • last_task_info: Dictionary caching inputs, initial conditions, trajectories and metadata.
RAISES DESCRIPTION
ImportError

If the external SSA package StochasticSimulationsNew is not available.
ValueError

If the network contains unknown parameters (num_unknown_params > 0) and the backend requires fully specified propensities (backend-dependent).
Notes
  • This method assumes the SSA helper RL4CRN.utils.stochastic.quick_measurement_SSA returns a summary dataframe with columns compatible with the parsing logic in this implementation.
  • last_task_info uses keys: 'inputs', 'initial conditions', 'time_horizon', 'trajectories', 'trajectories_std', 'outputs', 'outputs_std', 'has_diverged'.

transient_response(u_list, x0_list, time_horizon, LARGE_NUMBER=10000.0)

Compute deterministic transient responses (ODE integration).

This method integrates the ODE dynamics:

\[\dot{x}(t) = S\,a(x(t), u),\]

for each pair (u, x0) in the Cartesian product u_list × x0_list, where inputs u are treated as constant in time for each scenario.

Results are cached in last_task_info with type == 'transient response'. If the cache already contains a transient response, the stored results are returned immediately.

Numerical stability

If the integrator fails or the solution becomes unstable (detected by an event when max(|x|) exceeds LARGE_NUMBER or becomes non-finite), the remaining portion of the trajectory is filled with LARGE_NUMBER.

PARAMETER DESCRIPTION
u_list

List of constant input vectors, each array-like of shape (p,).
x0_list

List of initial condition vectors, each array-like of shape (n,).
time_horizon

1D array-like of time points of shape (T,) at which to evaluate/interpolate the solution.
LARGE_NUMBER

Threshold used both as an instability detection bound and as a fill value when integration fails.

DEFAULT: 10000.0

RETURNS DESCRIPTION

Tuple (time_horizon, x_list, y_list, last_task_info):

  • time_horizon: 1D numpy array of shape (T,).
  • x_list: List of full state trajectories, one per scenario, each of shape (n, T).
  • y_list: List of output trajectories, one per scenario, each of shape (q, T), extracted via output_idx.
  • last_task_info: Dictionary caching the results. Keys include 'inputs', 'initial conditions', 'time_horizon', 'trajectories', 'outputs'.
RAISES DESCRIPTION
ValueError

If num_unknown_params > 0.
ValueError

If solver is not one of 'LSODA' or 'CVODE'.
Notes
  • For 'LSODA', integration is performed using SciPy's solve_ivp.
  • For 'CVODE', the method delegates to solve_with_cvode and interpolates the solver output onto time_horizon.

transient_response_piecewise(u_nested_list, x0_list, nested_time_horizon, LARGE_NUMBER=10000.0, force=False)

Compute deterministic transient responses with piecewise-constant inputs.

This method generalizes transient_response to input sequences. Each element of u_nested_list is a sequence of input vectors [u_0, ..., u_{K-1}] applied consecutively over corresponding time segments nested_time_horizon = [t_0, ..., t_{K-1}].

For each scenario (u_sequence, x0) in the Cartesian product u_nested_list × x0_list, the method integrates the dynamics on each time segment in order, using the final state of segment k as the initial condition for segment k+1.

Instability handling

If any segment fails or terminates early due to instability, the remainder of the trajectory (for that scenario) is filled with LARGE_NUMBER.

Caching

Results are cached in last_task_info with type == 'transient response' and returned if a matching cache is detected (currently a heuristic check based on lengths).

PARAMETER DESCRIPTION
u_nested_list

List of input sequences. Each entry is a list of input vectors [u_0, ..., u_{K-1}], each vector array-like of shape (p,). The sequence length must match len(nested_time_horizon).
x0_list

List of initial conditions, each array-like of shape (n,).
nested_time_horizon

List of 1D time grids [t_0, ..., t_{K-1}]. Each t_k is a 1D array of times for segment k. The segments are concatenated (with offsets) into a single full_time_horizon stored in the cache and returned.
LARGE_NUMBER

Instability bound and fill value for failed segments.

DEFAULT: 10000.0

RETURNS DESCRIPTION

Tuple (full_time_horizon, x_list, y_list, last_task_info):
  • full_time_horizon: 1D array formed by concatenating shifted segment time grids, shape (T_total,).
  • x_list: List of full state trajectories, one per scenario, each of shape (n, T_total).
  • y_list: List of output trajectories, one per scenario, each of shape (q, T_total).
  • last_task_info: Cache dictionary storing inputs, initial conditions, time horizon, trajectories, and outputs.
RAISES DESCRIPTION
ValueError

If any input sequence length does not match the number of time segments.
ValueError

If num_unknown_params > 0.
ValueError

If solver is not one of 'LSODA' or 'CVODE'.

plot_logic_response(*, u_list=None, target_fn=None, title='Truth table', figsize=(2.2, 1.2), silent=False)

Plot a truth-table style view of the CRN's logic behavior.

This expects the CRN to have cached a deterministic multi-scenario simulation in self.last_task_info (e.g., produced by a reward function that calls transient_response(...) and stores its results).

PARAMETER DESCRIPTION
u_list

Optional list of input vectors used for the scenarios. If None, tries to read self.last_task_info['u_list'].

DEFAULT: None

target_fn

Optional callable mapping u (np.ndarray) -> {0,1} (or bool). If None, tries to read self.last_task_info['logic_fn'].

DEFAULT: None

title

Plot title.

DEFAULT: 'Truth table'

figsize

Matplotlib figure size.

DEFAULT: (2.2, 1.2)

silent

If True, does not call plt.show().

DEFAULT: False

RETURNS DESCRIPTION

(fig, ax) from plot_truth_table_transposed_nature.
RAISES DESCRIPTION
ValueError

If required cached information is missing.

plot_transient_response(fig=None, axes=None, alpha=0.1, plot_cfg=None)

Plot cached deterministic transient response trajectories.

This method visualizes trajectories stored by a call that cached last_task_info['type'] == 'transient response', typically produced by transient_response(...) or transient_response_piecewise(...).

Paper-ready styling can be controlled through plot_cfg, which may override default rcParams and common cosmetics (grid/spines/legend) while preserving the original method API.

PARAMETER DESCRIPTION
fig

Optional matplotlib Figure to draw into. If None and axes is None, a new Figure is created.

DEFAULT: None

axes

Optional axes list/array. If None and fig is None, new axes are created with one subplot per output.

DEFAULT: None

alpha

Line transparency for overlaid trajectories (default: 0.1).

DEFAULT: 0.1

plot_cfg

Optional dict controlling paper-ready styling. Supported keys: - rc: dict of matplotlib rcParams overrides - figsize: tuple (w,h) for figure size - constrained_layout: bool - tight_layout: bool - grid: bool - grid_kwargs: dict passed to ax.grid(...) - despine: bool (hide top/right spines) - spine_lw: float - alpha: override alpha passed in args - lw: line width override - title/xlabel/ylabel: override strings (title can be per-axis via None) - legend: "auto" | True | False - legend_kwargs: dict passed to ax.legend(...) - save: dict like {"path": "...pdf", "dpi": 600}

DEFAULT: None

RETURNS DESCRIPTION

Tuple[Figure, List[Axes]]: (fig, axes) where axes has length num_outputs.
RAISES DESCRIPTION
ValueError

If no cached deterministic transient response is present.

plot_transient_response_piecewise(fig=None, axes=None, alpha=1.0, input_color='red', shade_on=True, trj_color='green', pulse_lw=1.0, pulse_ls='--', gap=False, plot_cfg=None, is_main=True, normalize=False, y_label=None, legend_label=None)

Plot cached deterministic piecewise transient responses.

If is_main=False, this acts as a follower plotter: - plots only trajectories - uses the provided alpha - skips pulse overlays, shading, titles, and styling extras

If normalize=True, all trajectories are divided by a single global peak (max value across all cached conditions/scenarios/outputs in this object).

plot_phase_portrait(fig=None, axis=None, alpha=0.1, plot_cfg=None)

Plot a phase portrait from cached deterministic transient trajectories.

For two species, plots x1(t) vs x2(t). For three species, plots a 3D trajectory.

PARAMETER DESCRIPTION
fig

Optional matplotlib Figure.

DEFAULT: None

axis

Optional matplotlib Axis (2D or 3D).

DEFAULT: None

alpha

Line transparency for overlaid trajectories.

DEFAULT: 0.1

plot_cfg

Optional dict controlling paper-ready styling. Supported keys are the same as in plot_transient_response.

DEFAULT: None

RETURNS DESCRIPTION

Tuple[Figure, Axes]: (fig, axis).
RAISES DESCRIPTION
ValueError

If no cached deterministic transient response is present.
ValueError

If num_species is not 2 or 3.

plot_dose_response(fig=None, axes=None, alpha=0.5, plot_cfg=None)

Plot dose-response curves from cached simulation data.

Two cache modes are supported:

1) Dose-response cache: - last_task_info['type'] == 'dose response' - last_task_info['input doses'] : array-like - last_task_info['outputs'] : list of arrays shaped (q, n_doses) or similar

2) Transient-response-derived dose response: - last_task_info['type'] == 'transient response' - last_task_info['inputs'] : list of u vectors (assumes scalar dose is u[0]) - last_task_info['initial_conditions'] : list of x0 vectors - last_task_info['outputs'] : list of arrays shaped (q, T) for each (ic,u) scenario The method groups outputs by IC block and plots final-time output vs u[0].

Paper-ready styling can be controlled via plot_cfg.

PARAMETER DESCRIPTION
fig

Optional matplotlib Figure.

DEFAULT: None

axes

Optional axes list/array. If None, new axes are created (one per output).

DEFAULT: None

alpha

Line transparency for curves (default: 0.5).

DEFAULT: 0.5

plot_cfg

Optional dict controlling paper-ready styling (see plot_transient_response).

DEFAULT: None

RETURNS DESCRIPTION

Tuple[Figure, List[Axes]]: (fig, axes).
RAISES DESCRIPTION
ValueError

If neither dose-response nor transient-response cache is present.
KeyError

If required cache keys are missing.

plot_frequency_content(fig=None, axes=None, alpha=0.1, t0=0.0, plot_cfg=None)

Plot normalized one-sided FFT magnitude spectra from cached transients.

For each cached trajectory: 1) keep samples with time >= t0, 2) subtract mean to remove DC component, 3) compute one-sided FFT magnitudes, 4) normalize each spectrum by its maximum, 5) overlay spectra across scenarios.

PARAMETER DESCRIPTION
fig

Optional matplotlib Figure.

DEFAULT: None

axes

Optional axes list/array. If None, new axes are created (one per output).

DEFAULT: None

alpha

Line transparency for overlaid spectra.

DEFAULT: 0.1

t0

Time cutoff; only samples with time >= t0 are used.

DEFAULT: 0.0

plot_cfg

Optional dict controlling paper-ready styling (see plot_transient_response).

DEFAULT: None

RETURNS DESCRIPTION

Tuple[Figure, List[Axes]]: (fig, axes).
RAISES DESCRIPTION
ValueError

If no cached deterministic transient response is present.
ValueError

If insufficient points exist after t0.
ValueError

If inferred sampling interval is non-positive.

plot_SSA_transient_response(fig=None, axes=None, alpha=0.2, plot_cfg=None)

Plot cached SSA transient responses as mean ± std bands.

Expects stochastic cache produced by transient_response_SSA(...): - last_task_info['type'] == 'transient response SSA' - last_task_info['time_horizon'] : np.ndarray - last_task_info['outputs'] : List[np.ndarray] mean trajectories (q, T) - last_task_info['outputs_std'] : List[np.ndarray] std trajectories (q, T) - optionally last_task_info['inputs'] : List[u] for labeling

PARAMETER DESCRIPTION
fig

Optional matplotlib Figure.

DEFAULT: None

axes

Optional axes list/array. If None, new axes are created (one per output).

DEFAULT: None

alpha

Transparency for std shading (default: 0.2).

DEFAULT: 0.2

plot_cfg

Optional dict controlling paper-ready styling (see plot_transient_response). Note: plot_cfg["alpha"] (if set) overrides the line alpha, while the std-band alpha defaults to the method argument alpha unless plot_cfg includes "band_alpha".

Additional supported key: - band_alpha: float to override std-band transparency

DEFAULT: None

RETURNS DESCRIPTION

Tuple[Figure, List[Axes]]: (fig, axes).
RAISES DESCRIPTION
ValueError

If no SSA cache is present.
KeyError

If required cache keys are missing.

solve_with_cvode(x0, time_horizon, u, nonneg_idx, stop_fn)

Integrate the IOCRN ODE using CVODE and return a solve_ivp-like solution.

This is an internal helper used when self.solver == 'CVODE'. It wraps the sksundae.cvode.CVODE interface into an object that mimics SciPy's solve_ivp result structure (fields .t, .y, .status, .message, and .raw).

PARAMETER DESCRIPTION
x0

Initial condition vector of shape (n,).
time_horizon

1D array of requested times at which the solver should return values (similar to t_eval in SciPy).
u

Constant input vector of shape (p,) used by the rate function.
nonneg_idx

Indices of state variables constrained to be nonnegative. If provided, CVODE inequality constraints are applied.
stop_fn

Event-like function stop_fn(t, x) returning a scalar whose sign determines whether integration should stop (similar to SciPy event functions). Its attributes terminal and direction (if present) are forwarded to the CVODE events wrapper.
RETURNS DESCRIPTION

A lightweight solution object with:

  • t: 1D array of solver times
  • y: 2D array of shape (n, len(t))
  • status: integer status code
  • message: solver message
  • raw: the underlying CVODE solution object
RAISES DESCRIPTION
ImportError

If sksundae/CVODE is not installed.