ShapCRN

There is no learning without having to pose a question.
And a question requires doubt.
— Richard Feynman

Introduction

ShapCRN is a command-line application for studying SBML biochemical reaction-network models through simulation and controlled perturbations. It is designed to support both exploratory analysis and reproducible experiments: you can run model dynamics over time, inspect behavior near steady state, generate publication-friendly outputs (CSV and plots), and compare how system behavior changes when species or reactions are altered.

Beyond plain simulation, the project provides analysis pipelines to quantify influence and robustness at network level. In practice, this includes knockout/knockin workflows, Shapley-style importance assessment (with optional random or fixed perturbation scenarios), and Sobol-based global sensitivity analysis for selected targets. The file src/shapcrn/examples/usage_example.py is the CLI runner used by both module execution and the shapcrn console command.

Table of contents

• Overview
• Architecture and code map
• Main functionalities
• Additional functionalities
• Requirements
• Quickstart
• Commands
• Python API
• Output structure
• Examples
• Tips and troubleshooting
• Project layout

Overview

This project lets you load an SBML model and run one of several workflows:

• simulate: run time-course simulation with selectable integrators, optional steady-state mode, and CSV/plot exports.
• importance_assessment: estimate influence of species/reactions with knockout or knockin scenarios using a Shapley-style workflow.
• importance_assessment with perturbations: compare behavior under random or fixed perturbations of selected input species.
• sensitivity_analysis: compute Sobol indices for selected targets and optionally run convergence checks.
• knockout_species: create and save a modified SBML model where one species is disabled.
• knockout_reaction: create and save a modified SBML model where one reaction is disabled.
• knockin_species: create and save a modified SBML model where one species is reinforced/activated.
• knockin_reaction: create and save a modified SBML model where one reaction is reinforced/activated.

In a typical run, the flow is:

1. Parse command-line arguments (src/shapcrn/utils/utils.py).
2. Load and prepare SBML (src/shapcrn/utils/sbml/io.py + src/shapcrn/utils/sbml/reactions.py).
3. Dispatch to a pipeline (src/shapcrn/pipelines/...).
4. Run simulations/analysis (src/shapcrn/utils/simulation.py, src/shapcrn/utils/sensitivity.py).
5. Save artifacts (CSV, plots, reports, edited SBML).
6. Return logs and outputs under the selected output folder.

SBML compatibility and preparation

Models keep their original SBML level/version. Core transformations are tested with Level 2 Versions 1 and 4 and Level 3 Versions 1 and 2. Loading, preparation, and saving reject libSBML errors with validation codes and details. Warnings, such as unspecified units, remain nonfatal and are written when a log path is provided. Invalid output is rejected before an existing output file is opened.

Reaction splitting retains the project's kinetic-inference policy, independently of the reversible flag. Expressions containing division retain the historical Michaelis–Menten exclusion. Supported candidates are a binary difference such as kf*A - kr*B, optionally multiplied by common factors, or an equivalent function body. AST decomposition preserves the original expression: cell*kf*A - kr*B becomes forward cell*kf*A and reverse kr*B; it does not add a compartment factor to the reverse rate. Function arguments retain their order and scope, and original function definitions remain available to other callers.

Splitting preserves local parameter units, modifiers and fixed stoichiometry. Generated L3V1 reactions have fast=false; L3V2 reactions omit fast. The constant attribute on a species reference describes fixed stoichiometry, independently of the referenced species' constant flag. The low-level create_sbml_reaction_LMA helper accepts (species_id, stoichiometry, constant) triples for reactants/products and species IDs for modifiers. Creation and split helpers return objects already attached to the model; callers must not add them again. Legacy splitting retains its _forward/_reverse suffixes.

An unsupported inferred candidate raises an error identifying the reaction. Dynamic or assignment-defined stoichiometry, dependencies on replaced reaction or species-reference IDs, fast=true transformations, and populated SBML packages are rejected before modifying the input. Empty legacy Layout plugins attached by libSBML to ordinary Level 2 files are allowed. A failed transformation leaves the input model unchanged. SBML permits missing kinetics, but preparation and rate-based analysis require them. This support is not a general-purpose splitter for every mathematically valid SBML kinetic law.

Sampling, sensitivity inputs, peak detection and knock-in values use the units of the SBML species symbol: amounts when hasOnlySubstanceUnits=true, and concentrations otherwise. This is independent of whether the XML stores initialAmount or initialConcentration; conversion uses compartment size. For example, initialAmount=12 in a compartment of size 4 corresponds to a sample value of 3 when hasOnlySubstanceUnits=false. Unresolved initialization or an unavailable conversion is reported instead of supplying a numeric default. Existing callers that supplied stored amounts for concentration-valued symbols must now supply concentrations. Explicit simulation output selections remain under the caller's control.

Reaction knock-in creates fixed species copies and replaces all applicable kinetic-law occurrences while respecting local parameter shadowing. Repeated references to one species share a copy and require identical input values. Species knock-in rejects conflicting rules, events and initial assignments; species knockout instead sets its existing updates to zero. Knockout retains valid participants and modifier declarations when removing product references. The deactivation helper creates a rapid ordinary sink, with concentration/amount and conversion-factor handling; it does not assert SBML's instantaneous fast approximation.

Specification references: SBML L3V1 Core, §3.3 (identifiers), §4.3 (functions), §4.6 (species quantities), and §4.11 (reactions, stoichiometry, modifiers and local parameters).

Regression tests can be run with:

conda run -n shapcrn-dev python -m pytest -q

Architecture and code map

The codebase follows a layered structure:

• Example runner layer: src/shapcrn/examples/usage_example.py demonstrates how to wire commands to pipelines. It is also the CLI entry implementation used by python -m ... and by the shapcrn console script.
• Pipeline layer: src/shapcrn/pipelines/* contains command-oriented orchestration. Each pipeline parses command-specific arguments, coordinates utilities, and writes outputs.
• Utility layer: src/shapcrn/utils/* contains reusable logic split by domain:
  • src/shapcrn/utils/sbml/: SBML I/O, reaction preprocessing, and knock operations.
  • src/shapcrn/utils/simulation.py: RoadRunner setup, simulation, perturbation sampling and aggregation helpers.
  • src/shapcrn/utils/sensitivity.py: Sobol setup/execution, convergence checks, and statistics.
  • src/shapcrn/utils/plot.py: static and interactive plotting utilities.
  • src/shapcrn/utils/graph.py: model-to-network conversion and graph rendering.
  • src/shapcrn/utils/utils.py: CLI parser construction, normalization helpers, logging, and shared helpers.

Design intent:

• Keep pipelines thin and scenario-focused.
• Keep model/math/plot logic reusable in utils.
• Keep output organization consistent across commands (images/, csv/, reports/, dot/).

Main functionalities

This section describes the core capabilities and maps them to the main functions in the codebase.

1) Simulate

Main execution path:

• Example runner: src/shapcrn/examples/usage_example.py (command simulate)
• Model loading and normalization: src/shapcrn/utils/sbml/io.py::load_and_prepare_model
• RoadRunner setup: src/shapcrn/utils/simulation.py::load_roadrunner_model
• Simulation engine: src/shapcrn/utils/simulation.py::simulate
• Steady-state mode (if enabled): src/shapcrn/utils/simulation.py::simulate_with_steady_state
• Output plots: src/shapcrn/utils/plot.py::plot_results and plot_results_interactive

How it works:

• load_roadrunner_model(...) converts the SBML model to a RoadRunner instance and configures integrator/tolerances.
• simulate(...) supports two modes:
  • Standard mode: one run from start_time to end_time.
  • Steady-state mode: adaptive block simulation until variation is below threshold.
• simulate_with_steady_state(...) compares the last point of consecutive blocks and tracks relative/absolute variation per monitored species. When all monitored species stay below threshold for consecutive checks, it flags steady state.
• Results are exported to CSV and plotted (static PNG or interactive HTML).

2) Species/Reactions Knockout

Main execution path:

• Species KO pipeline: src/shapcrn/pipelines/knockout/knockout_species.py::knockout_species
• Reaction KO pipeline: src/shapcrn/pipelines/knockout/knockout_reaction.py::knockout_reaction
• Core KO logic: src/shapcrn/utils/sbml/knock.py::knockout_species and knockout_reaction

How species knockout works (knockout_species in knock.py):

• Forces assignment/initial rules for the target species to 0 when present.
• Scans events and initial assignments and sets the target species update math to 0.
• If the species is a reactant in a reaction, that reaction is marked for knockout.
• If the species is only a product, the product entry is removed from that reaction; if no products remain, the reaction is also knocked out.
• Finalizes by setting target species initial concentration to 0.0 and boundary condition True (fixed species).

How reaction knockout works (knockout_reaction in knock.py):

• Looks up the target reaction.
• Replaces the kinetic law AST with constant 0, disabling flux while preserving the reaction object in the SBML structure.

3) Species/Reactions Knockin

Main execution path:

• Species KI pipeline: src/shapcrn/pipelines/knockin/knockin_species.py::knockin_species
• Reaction KI pipeline: src/shapcrn/pipelines/knockin/knockin_reaction.py::knockin_reaction
• Core KI logic: src/shapcrn/utils/sbml/knock.py::knockin_species and knockin_reaction

How species knockin works:

• get_species_peak_value(...) runs a short simulation (end_time=60) and uses the maximum value in the species symbol's amount/concentration units as the knock-in value.
• knockin_species(...) then sets that value according to hasOnlySubstanceUnits and marks the species as fixed (boundaryCondition=True, constant=True).

How reaction knockin works:

• get_reactants_peak_values(...) collects max simulated values for each reactant of the target reaction.
• knockin_reaction(...) creates constant reactant copies, one per original reactant.
• The target reaction is cloned, its reactants are replaced with these new constant species, and the kinetic law expression is rewritten to reference species.
• The original reaction is removed and the modified cloned reaction is added back to the model.

Additional functionalities

Importance assessment

Main path:

• Pipeline: src/shapcrn/pipelines/importance.py::importance_assessment
• Model setup: model_preparation(...)
• Sampling setup: generate_samples(...)
• Baseline simulation: simulate_original_model(...)
• KO/KI simulation batch: simulate_knocked_data(...)
• Payoff/Shapley: run_shap_analysis(...)
• Perturbation diagnostics: assess_perturbation_importance(...) and report generation

What it produces:

• Variation matrices (log-ratio based)
• Shapley matrices (raw and normalized for plotting)
• Heatmaps and optional text reports
• Optional fixed-vs-random perturbation comparison outputs

Sensitivity analysis

Main path:

• Pipeline: src/shapcrn/pipelines/sensitivity_analysis.py::sensitivity_analysis
• Problem specification: src/shapcrn/utils/sensitivity.py::get_problem_parameters
• Sobol sampling/analysis: SALib (sobol.sample, sobol.analyze)
• Batch simulation backend: src/shapcrn/utils/sensitivity.py::run_simulation_with_params
• Optional convergence workflow: run_convergence_analysis(...) + convergence plots

What it produces:

• Sobol indices per target species
• Optional convergence diagnostics and plots
• Fixed-vs-sampled perturbation comparison CSV (when requested)

Network generation

Main path:

• Pipeline: src/shapcrn/pipelines/network.py::create_model_network
• Graph extraction: src/shapcrn/utils/graph.py::get_network_from_sbml
• Rendering: src/shapcrn/utils/graph.py::plot_network

What it produces:

• Network image (PNG)
• DOT graph source (.gv) for external graph tooling

Requirements

ShapCRN supports Python 3.10, 3.11, and 3.12.

Core runtime dependencies:

• numpy
• pandas
• scipy
• scikit-learn
• python-libsbml
• libroadrunner
• SALib
• matplotlib
• seaborn
• plotly
• networkx

Install the released package:

python -m pip install shapcrn

Network rendering uses PyGraphviz and the Graphviz system package. Install the optional Python dependency with:

python -m pip install "shapcrn[network]"

For development from a source checkout:

python -m pip install -e ".[dev,network]"

Quickstart

After installation, check the CLI:

shapcrn -h

You can also run the module entrypoint directly:

python -m shapcrn.examples.usage_example -h

The following command assumes a source checkout. For an installed package, replace the model path with your own SBML file:

shapcrn simulate models/KnockinModelV2.xml -t 120 -o results

Inspect command-specific options:

python -m shapcrn.examples.usage_example simulate -h
python -m shapcrn.examples.usage_example importance_assessment -h
python -m shapcrn.examples.usage_example sensitivity_analysis -h

Commands

The shapcrn console command is canonical. The historical shapcrn.examples.usage_example module remains as a compatibility wrapper.

General form:

python -m shapcrn.examples.usage_example <command> [options]

Equivalent console-script form after install:

shapcrn <command> [options]

Available commands:

• simulate Time-course model simulation with optional steady-state termination and plotting.
• importance_assessment Shapley-style influence analysis with knockout/knockin scenarios, optionally with perturbations.
• sensitivity_analysis Global Sobol sensitivity analysis for selected species/targets.
• knockout_species Build and save a model where one species is disabled.
• knockout_reaction Build and save a model where one reaction is disabled.
• knockin_species Build and save a model where one species is fixed to a computed reinforced value.
• knockin_reaction Build and save a model where one reaction is reinforced via reactant replacement strategy.
• create_network Build and save a reaction-network graph from the SBML model.

Python API

The high-level API returns typed result objects and does not write files unless an output directory or output path is explicitly supplied:

from shapcrn import analyze_sensitivity, assess_importance, simulate_model

simulation = simulate_model("model.xml", end_time=120, points=500)
print(simulation.data.tail())

importance = assess_importance(
    "model.xml",
    operation="knockout",
    input_species=["S1", "S2"],
    use_perturbations=True,
    seed=42,
)

sensitivity = analyze_sensitivity(
    "model.xml",
    input_species=["S1", "S2"],
    base_samples=1024,
    seed=42,
)

In-memory SBML helpers and low-level simulation functions remain available:

from shapcrn import load_model, load_roadrunner_model, simulate
from shapcrn import knockout_species, save_sbml_model

Output structure

By default, outputs are written under ./results in a model-specific folder:

<output>/<model_name>/
├── csv/
├── dot/
├── images/
└── reports/

Simulation writes simulation.csv plus a PNG or HTML plot. Sensitivity writes sobol_indices.csv, sobol_interactions.csv, and a text report; convergence and fixed-perturbation artifacts are added only when requested.

Examples

1) Simulate a model (static plot + CSV)

python -m shapcrn.examples.usage_example simulate models/KnockinModelV2.xml \
  -t 120 \
  -i cvode \
  -o results

2) Simulate until steady state (interactive HTML plot)

python -m shapcrn.examples.usage_example simulate models/KnockinModelV2.xml \
  --steady-state \
  --max-time 2000 \
  --sim-step 10 \
  --threshold 1e-7 \
  --interactive \
  -o results

3) Importance assessment (knockout, no perturbations)

python -m shapcrn.examples.usage_example importance_assessment models/KnockinModelV2.xml \
  --operation knockout \
  --payoff-function last \
  -t 120 \
  -o results

4) Importance assessment with random perturbations

python -m shapcrn.examples.usage_example importance_assessment models/KnockinModelV2.xml \
  --operation knockout \
  --input-species S1 S2 \
  --use-perturbations \
  --num-samples 10 \
  --max-combinations 2000 \
  --variation 20 \
  --payoff-function max \
  -t 120 \
  -o results

5) Importance assessment with fixed perturbations

python -m shapcrn.examples.usage_example importance_assessment models/KnockinModelV2.xml \
  --operation knockin \
  --input-species S1 S2 \
  --use-perturbations \
  --use-fixed-perturbations \
  --fixed-perturbations -20 20 \
  --payoff-function min \
  -t 120 \
  -o results

6) Sensitivity analysis (Sobol)

python -m shapcrn.examples.usage_example sensitivity_analysis models/KnockinModelV2.xml \
  --input-species S1 S2 \
  --base-samples 1024 \
  --perturbation-range 20 \
  --seed 42 \
  -o results

7) Sensitivity analysis with convergence check

python -m shapcrn.examples.usage_example sensitivity_analysis models/KnockinModelV2.xml \
  --input-species S1 S2 \
  --check-convergence \
  -o results

8) Knock out one species and save edited model

python -m shapcrn.examples.usage_example knockout_species models/KnockinModelV2.xml S1 \
  --model-dir models \
  -o results

9) Knock out one reaction and save edited model

python -m shapcrn.examples.usage_example knockout_reaction models/KnockinModelV2.xml R1_MassAction_Explicit \
  --model-dir models \
  -o results

10) Knock in one species and save edited model

python -m shapcrn.examples.usage_example knockin_species models/KnockinModelV2.xml S1 \
  --model-dir models \
  -o results

11) Knock in one reaction and save edited model

python -m shapcrn.examples.usage_example knockin_reaction models/KnockinModelV2.xml R1_MassAction_Explicit \
  --model-dir models \
  -o results

Tips and troubleshooting

• models/KnockinModelV2.xml works with the IDs used in these examples (S1, S2, R1_MassAction_Explicit).
• The files under models/ are repository examples and are intentionally not included in PyPI artifacts.
• Add logging to any command with -l <log_file_path>.
• For large perturbation spaces, use --max-combinations to cap Cartesian-product runs and avoid RAM saturation.
• If simulation fails immediately, first run simulate -h and verify required options and valid integrator/model values.
• If you see ModuleNotFoundError: No module named 'shapcrn', run python -m pip install -e . from repo root or run commands with PYTHONPATH=src.

Project layout

• pyproject.toml: package metadata and console-script definition
• src/shapcrn/: Python package root
• src/shapcrn/api.py: stable programmatic API and result types
• src/shapcrn/cli.py: console command adapter
• src/shapcrn/examples/usage_example.py: compatibility entry point
• models/: sample SBML models
• results/: default output directory

License

ShapCRN is distributed under the MIT License. See LICENSE.

Acknowledgments

To Aurora
Thank you for your constant support, encouragement, and guidance.
This work wouldn't be what it is without you.

---