libsbml-network-modeling

Overview

libSBML is the reference library for reading, writing, creating, and validating SBML (Systems Biology Markup Language) models. SBML is the community standard for encoding biochemical reaction networks — ODE models, signaling cascades, and genome-scale metabolic models all use it. The Python API (python-libsbml) exposes a full object model covering compartments, species, reactions, kinetic laws, rules, constraints, and every SBML extension. Models saved as SBML .xml files are interoperable with COPASI, Tellurium, RoadRunner, COBRApy, and BioModels Database.

When to Use

• Building a new ODE-based biochemical model (enzyme kinetics, signaling pathway) from scratch in SBML format for simulation in COPASI or Tellurium
• Reading and programmatically modifying an existing BioModels Database model — changing kinetic parameters, adding species, or patching reaction stoichiometry
• Validating an SBML file against the specification before submitting to BioModels or sharing with collaborators
• Converting SBML models between Level 1/2/3 for compatibility with older simulation tools
• Constructing genome-scale metabolic models with flux bounds and objective functions via the FBC (Flux Balance Constraints) extension for use with COBRApy
• Parsing an SBML model to extract the stoichiometry matrix, species list, or reaction network as NumPy/pandas data structures for custom analysis
• Use cobrapy-metabolic-modeling instead when you need to run FBA, FVA, or gene knockouts on an already-built metabolic model — libSBML is for constructing and editing the SBML file itself
• Use tellurium directly when you want an integrated Python environment for both SBML authoring (Antimony syntax) and ODE simulation without low-level XML manipulation

Prerequisites

• Python packages: python-libsbml, numpy, pandas (optional, for matrix extraction)
• Optional packages: cobra (COBRApy, for FBA after SBML load), tellurium (for SBML↔Antimony conversion and simulation)
• Data requirements: SBML files (.xml), or built from scratch in Python; BioModels Database SBML files are freely available at https://www.ebi.ac.uk/biomodels/

pip install python-libsbml numpy pandas
# Optional simulation/FBA integrations:
pip install cobra tellurium

Quick Start

Load an SBML file, inspect its content, and modify a parameter value:

import libsbml

# Read an SBML model file
reader = libsbml.SBMLReader()
doc = reader.readSBMLFromFile("BIOMD0000000012.xml")

# Check for errors
if doc.getNumErrors() > 0:
    doc.printErrors()

model = doc.getModel()
print(f"Model: {model.getId()}")
print(f"  Compartments: {model.getNumCompartments()}")
print(f"  Species:      {model.getNumSpecies()}")
print(f"  Reactions:    {model.getNumReactions()}")

# Modify a global parameter
param = model.getParameter("Km")
if param:
    old_val = param.getValue()
    param.setValue(0.05)
    print(f"Updated Km: {old_val} → {param.getValue()}")

# Write modified model back to file
writer = libsbml.SBMLWriter()
writer.writeSBMLToFile(doc, "BIOMD0000000012_modified.xml")
print("Saved modified model.")

Core API

Module 1: Reading and Validating SBML

Load SBML files from disk or strings, check parse errors, and run full SBML spec validation.

import libsbml

# Read from file
reader = libsbml.SBMLReader()
doc = reader.readSBMLFromFile("model.xml")

# Check for fatal parse errors
n_errors = doc.getNumErrors()
print(f"Parse errors: {n_errors}")
for i in range(n_errors):
    err = doc.getError(i)
    severity = err.getSeverityAsString()
    print(f"  [{severity}] line {err.getLine()}: {err.getMessage()}")

# Check the SBML Level and Version
print(f"SBML Level {doc.getLevel()} Version {doc.getVersion()}")

# Read from in-memory XML string
xml_string = open("model.xml").read()
doc2 = reader.readSBMLFromString(xml_string)
model = doc2.getModel()
print(f"Model id: {model.getId()}, name: {model.getName()}")

import libsbml

# Full consistency / validation check (more thorough than parse error check)
doc = libsbml.readSBMLFromFile("model.xml")

# Enable all consistency checks
doc.setConsistencyChecks(libsbml.LIBSBML_CAT_GENERAL_CONSISTENCY, True)
doc.setConsistencyChecks(libsbml.LIBSBML_CAT_IDENTIFIER_CONSISTENCY, True)
doc.setConsistencyChecks(libsbml.LIBSBML_CAT_UNITS_CONSISTENCY, True)
doc.setConsistencyChecks(libsbml.LIBSBML_CAT_MATHML_CONSISTENCY, True)
doc.setConsistencyChecks(libsbml.LIBSBML_CAT_SBO_CONSISTENCY, True)
doc.setConsistencyChecks(libsbml.LIBSBML_CAT_OVERDETERMINED_MODEL, True)
doc.setConsistencyChecks(libsbml.LIBSBML_CAT_MODELING_PRACTICE, True)

n_errors = doc.checkConsistency()
print(f"Consistency check: {n_errors} issue(s)")
for i in range(n_errors):
    err = doc.getError(i)
    print(f"  [{err.getSeverityAsString()}] {err.getShortMessage()}: {err.getMessage()[:120]}")

Module 2: Creating Models from Scratch

Build a complete SBML document by adding compartments, species, and reactions programmatically.

import libsbml

# Create a new SBML Level 3 Version 2 document
doc = libsbml.SBMLDocument(3, 2)
model = doc.createModel()
model.setId("simple_enzymatic_model")
model.setName("Simple Enzymatic Reaction Model")
model.setTimeUnits("second")
model.setSubstanceUnits("mole")
model.setVolumeUnits("litre")
model.setExtentUnits("mole")

# Add a compartment (cytoplasm)
comp = model.createCompartment()
comp.setId("cytoplasm")
comp.setName("Cytoplasm")
comp.setConstant(True)
comp.setSize(1.0)          # 1 litre
comp.setSpatialDimensions(3)
comp.setUnits("litre")

# Add species: substrate S, enzyme E, complex ES, product P
species_data = [
    ("S",  "Substrate",           0.01, True),   # (id, name, initialConc, boundaryCondition)
    ("E",  "Enzyme",              0.001, False),
    ("ES", "Enzyme-Substrate",    0.0,  False),
    ("P",  "Product",             0.0,  True),
]
for sp_id, sp_name, init_conc, boundary in species_data:
    sp = model.createSpecies()
    sp.setId(sp_id)
    sp.setName(sp_name)
    sp.setCompartment("cytoplasm")
    sp.setInitialConcentration(init_conc)
    sp.setBoundaryCondition(boundary)
    sp.setHasOnlySubstanceUnits(False)
    sp.setConstant(False)
    print(f"Added species: {sp_id} (init={init_conc} M, boundary={boundary})")

print(f"Model has {model.getNumSpecies()} species and {model.getNumCompartments()} compartment(s)")

Module 3: Editing Reactions and Kinetic Laws

Add reactions with stoichiometry and MathML kinetic law formulas.

import libsbml

# Continuing from Module 2: add Michaelis-Menten kinetics reactions
# Forward: S + E -> ES (association)
# Reverse: ES -> S + E (dissociation)
# Catalytic: ES -> P + E (product release)

# First, add kinetic parameters as global parameters
params = [
    ("kf", 1e6,  "litre per mole per second"),   # forward rate constant
    ("kr", 1e-3, "per second"),                   # reverse rate constant
    ("kcat", 0.1, "per second"),                  # catalytic rate constant
]
for p_id, p_val, p_units in params:
    param = model.createParameter()
    param.setId(p_id)
    param.setValue(p_val)
    param.setConstant(True)
    # Units are for documentation — libSBML stores them as unit definitions
    print(f"Added parameter: {p_id} = {p_val}")

def add_reaction(model, rxn_id, rxn_name, reactants, products, formula):
    """Helper: create a reaction with MathML kinetic law."""
    rxn = model.createReaction()
    rxn.setId(rxn_id)
    rxn.setName(rxn_name)
    rxn.setReversible(False)
    for sp_id, stoich in reactants:
        sr = rxn.createReactant()
        sr.setSpecies(sp_id)
        sr.setStoichiometry(stoich)
        sr.setConstant(True)
    for sp_id, stoich in products:
        sr = rxn.createProduct()
        sr.setSpecies(sp_id)
        sr.setStoichiometry(stoich)
        sr.setConstant(True)
    kl = rxn.createKineticLaw()
    math_ast = libsbml.parseL3Formula(formula)
    if math_ast is None:
        raise ValueError(f"Could not parse formula: {formula}")
    kl.setMath(math_ast)
    return rxn

add_reaction(model, "v1", "Association",  [("S",1),("E",1)], [("ES",1)], "kf * S * E * cytoplasm")
add_reaction(model, "v2", "Dissociation", [("ES",1)], [("S",1),("E",1)], "kr * ES * cytoplasm")
add_reaction(model, "v3", "Catalysis",    [("ES",1)], [("P",1),("E",1)], "kcat * ES * cytoplasm")

print(f"Model has {model.getNumReactions()} reactions")
writer = libsbml.SBMLWriter()
writer.writeSBMLToFile(doc, "michaelis_menten.xml")
print("Saved michaelis_menten.xml")

Module 4: Species and Compartments

Inspect and modify species properties — initial amounts vs concentrations, boundary conditions, compartment volumes.

import libsbml

doc = libsbml.readSBMLFromFile("model.xml")
model = doc.getModel()

# Iterate species and print their properties
print(f"{'ID':<12} {'Compartment':<15} {'InitConc':>10} {'InitAmt':>10} {'Boundary':>10} {'Constant':>10}")
print("-" * 70)
for i in range(model.getNumSpecies()):
    sp = model.getSpecies(i)
    init_conc = sp.getInitialConcentration() if sp.isSetInitialConcentration() else "—"
    init_amt  = sp.getInitialAmount()        if sp.isSetInitialAmount()        else "—"
    print(f"{sp.getId():<12} {sp.getCompartment():<15} {str(init_conc):>10} {str(init_amt):>10} "
          f"{str(sp.getBoundaryCondition()):>10} {str(sp.getConstant()):>10}")

# Modify compartment volume (e.g. scale to a smaller cell)
comp = model.getCompartment("cytoplasm")
if comp:
    old_size = comp.getSize()
    comp.setSize(old_size * 0.1)
    print(f"\nCytoplasm volume: {old_size} → {comp.getSize()} litre")

# Set a species initial concentration by ID
sp = model.getSpecies("S")
if sp:
    sp.setInitialConcentration(0.005)
    print(f"Updated [S] initial concentration to {sp.getInitialConcentration()} M")

Module 5: Rules and Constraints

Add assignment rules, rate rules, and algebraic rules to model derived quantities or conserved relationships.

import libsbml

doc = libsbml.readSBMLFromFile("model.xml")
model = doc.getModel()

# AssignmentRule: computes a variable algebraically at every time step
# Example: total enzyme E_total = E + ES (conservation relationship, monitoring only)
ar = model.createAssignmentRule()
ar.setVariable("E_total")   # must be an existing parameter or species id
# First add the target as a parameter if needed
if model.getParameter("E_total") is None:
    p = model.createParameter()
    p.setId("E_total")
    p.setConstant(False)     # MUST be False for assignment rule targets
    p.setValue(0.0)
math_ast = libsbml.parseL3Formula("E + ES")
ar.setMath(math_ast)
print(f"Added AssignmentRule: E_total = E + ES")

# RateRule: specifies dX/dt directly, bypassing reaction-based ODE generation
# Useful for custom non-mass-action dynamics
rr = model.createRateRule()
rr.setVariable("P")   # species P already has boundary=True to allow rate rules
math_ast2 = libsbml.parseL3Formula("kcat * ES * cytoplasm")
rr.setMath(math_ast2)
print(f"Added RateRule: dP/dt = kcat * ES * cytoplasm")

# Constraint: model invariant that simulators should monitor (not enforced computationally)
constraint = model.createConstraint()
math_ast3 = libsbml.parseL3Formula("S >= 0")
constraint.setMath(math_ast3)
msg = libsbml.XMLNode.convertStringToXMLNode("<message><p>Substrate cannot be negative</p></message>")
constraint.setMessage(msg)
print(f"Added Constraint: S >= 0")

print(f"Model rules: {model.getNumRules()}, constraints: {model.getNumConstraints()}")

Module 6: FBC Extension (Flux Balance Constraints)

Use the SBML FBC package to encode genome-scale metabolic models with flux bounds and an objective function for use with COBRApy or other FBA solvers.

import libsbml

# Build a minimal FBC-enabled model (3-reaction toy network)
doc = libsbml.SBMLDocument(3, 2)
# Enable FBC package (required)
doc.enablePackage(libsbml.FbcExtension.getXmlnsL3V1V2(), "fbc", True)
doc.setPackageRequired("fbc", False)

model = doc.createModel()
model.setId("toy_fba_model")
fbc_plugin = model.getPlugin("fbc")
fbc_plugin.setStrict(True)

# Add a compartment and species
comp = model.createCompartment()
comp.setId("c")
comp.setConstant(True)
comp.setSize(1.0)

for sp_id in ["A", "B", "C"]:
    sp = model.createSpecies()
    sp.setId(sp_id)
    sp.setCompartment("c")
    sp.setInitialAmount(0.0)
    sp.setBoundaryCondition(False)
    sp.setConstant(False)
    sp.setHasOnlySubstanceUnits(True)
    sp_fbc = sp.getPlugin("fbc")
    sp_fbc.setChemicalFormula("")

# Add flux bound parameters
bounds = {"lb_0": 0.0, "lb_neg1000": -1000.0, "ub_1000": 1000.0}
for b_id, b_val in bounds.items():
    p = model.createParameter()
    p.setId(b_id)
    p.setValue(b_val)
    p.setConstant(True)

def add_fbc_reaction(model, rxn_id, reactants, products, lb_id, ub_id):
    rxn = model.createReaction()
    rxn.setId(rxn_id)
    rxn.setReversible(lb_id == "lb_neg1000")
    rxn.setFast(False)
    for sp_id, stoich in reactants:
        sr = rxn.createReactant(); sr.setSpecies(sp_id); sr.setStoichiometry(stoich); sr.setConstant(True)
    for sp_id, stoich in products:
        sr = rxn.createProduct();  sr.setSpecies(sp_id); sr.setStoichiometry(stoich); sr.setConstant(True)
    rxn_fbc = rxn.getPlugin("fbc")
    rxn_fbc.setLowerFluxBound(lb_id)
    rxn_fbc.setUpperFluxBound(ub_id)
    return rxn

add_fbc_reaction(model, "r1", [("A", 1)], [("B", 1)], "lb_0", "ub_1000")
add_fbc_reaction(model, "r2", [("B", 1)], [("C", 1)], "lb_0", "ub_1000")
add_fbc_reaction(model, "r3", [("A", 1)], [],          "lb_0", "ub_1000")  # exchange

# Add objective function: maximize r2 flux
obj = fbc_plugin.createObjective()
obj.setId("maximize_r2")
obj.setType("maximize")
fbc_plugin.setActiveObjectiveId("maximize_r2")
flux_obj = obj.createFluxObjective()
flux_obj.setReaction("r2")
flux_obj.setCoefficient(1.0)

writer = libsbml.SBMLWriter()
writer.writeSBMLToFile(doc, "toy_fba.xml")
print(f"Saved toy_fba.xml (Level {doc.getLevel()} Version {doc.getVersion()}, FBC enabled)")
print(f"Reactions: {model.getNumReactions()}, Objective: maximize r2")

Module 7: Exporting and Level Conversion

Write models to file or string, convert between SBML levels, and export to Antimony notation via Tellurium.

import libsbml

doc = libsbml.readSBMLFromFile("michaelis_menten.xml")
model = doc.getModel()
print(f"Loaded: Level {doc.getLevel()}, Version {doc.getVersion()}")

# Write to XML string (useful for in-memory transmission)
writer = libsbml.SBMLWriter()
xml_string = writer.writeSBMLToString(doc)
print(f"XML string length: {len(xml_string)} characters")

# Convert Level 3 → Level 2 (for compatibility with older tools)
# SBMLDocument.setLevelAndVersion handles conversion automatically
props = libsbml.ConversionProperties()
props.addOption("setLevelAndVersion", True, "Convert level and version")
props.addOption("targetLevel",   2)
props.addOption("targetVersion", 4)
status = doc.convert(props)
if status == libsbml.LIBSBML_OPERATION_SUCCESS:
    writer.writeSBMLToFile(doc, "michaelis_menten_L2V4.xml")
    print(f"Converted to L2V4 → michaelis_menten_L2V4.xml")
else:
    print(f"Conversion failed with code: {status}")

# Export SBML to Antimony (human-readable) via Tellurium (optional)
try:
    import tellurium as te
    antimony_str = te.sbmlToAntimony(open("michaelis_menten.xml").read())
    print("Antimony notation:")
    print(antimony_str[:600])
    with open("michaelis_menten.ant", "w") as f:
        f.write(antimony_str)
    print("Saved michaelis_menten.ant")
except ImportError:
    print("tellurium not installed — skipping Antimony export")

Key Concepts

SBMLDocument, Model, and the Plugin Architecture

Every libSBML session starts with an SBMLDocument that owns exactly one Model. Extension packages (FBC, qual, layout, groups, distrib) are accessed as plugins retrieved via object.getPlugin("fbc"). Plugins are only available after enabling the package on the document with doc.enablePackage(...). Calling getPlugin on a document that has not enabled the package returns None.

import libsbml

doc = libsbml.readSBMLFromFile("iJO1366.xml")
model = doc.getModel()

# Check which packages are active
for i in range(doc.getNumPlugins()):
    pkg = doc.getPlugin(i)
    print(f"Package: {pkg.getPackageName()} (level {pkg.getLevel()})")

# Access FBC plugin
fbc_plugin = model.getPlugin("fbc")
if fbc_plugin:
    print(f"FBC strict mode: {fbc_plugin.getStrict()}")
    print(f"Objectives: {fbc_plugin.getNumObjectives()}")

MathML Formulas and the AST

Kinetic laws in SBML are stored as MathML. libSBML parses formula strings to an Abstract Syntax Tree (AST) using libsbml.parseL3Formula(string) and converts AST back to a string with libsbml.formulaToL3String(ast). Always check that parseL3Formula returns non-None before assigning to a kinetic law — a None return means parsing failed silently.

import libsbml

# Parse and inspect a kinetic formula
formula = "Vmax * S / (Km + S) * cytoplasm"
ast = libsbml.parseL3Formula(formula)
if ast is None:
    print("ERROR: formula could not be parsed")
else:
    print(f"Parsed formula: {libsbml.formulaToL3String(ast)}")
    print(f"AST root type: {ast.getType()}")  # e.g., AST_TIMES

# Retrieve a kinetic law formula from an existing reaction
doc = libsbml.readSBMLFromFile("model.xml")
model = doc.getModel()
rxn = model.getReaction(0)
if rxn and rxn.isSetKineticLaw():
    kl = rxn.getKineticLaw()
    formula_str = libsbml.formulaToL3String(kl.getMath())
    print(f"Reaction '{rxn.getId()}' kinetic law: {formula_str}")

Common Workflows

Workflow 1: Load BioModels Model, Modify Parameters, and Simulate

Goal: Download a BioModels model, adjust kinetic parameters, and run an ODE simulation with Tellurium/RoadRunner.

import libsbml
import urllib.request

# 1. Download SBML from BioModels Database (BIOMD0000000012 = Tyson 1991 cell cycle)
url = "https://www.ebi.ac.uk/biomodels/model/download/BIOMD0000000012?filename=BIOMD0000000012_url.xml"
urllib.request.urlretrieve(url, "BIOMD0000000012.xml")
print("Downloaded BIOMD0000000012.xml")

# 2. Load and inspect the model
doc = libsbml.readSBMLFromFile("BIOMD0000000012.xml")
model = doc.getModel()
print(f"Model: {model.getId()} | Level {doc.getLevel()} Version {doc.getVersion()}")
print(f"Species: {model.getNumSpecies()}, Reactions: {model.getNumReactions()}")

# 3. Print all global parameters and their values
print("\nGlobal parameters:")
for i in range(model.getNumParameters()):
    p = model.getParameter(i)
    print(f"  {p.getId():<20} = {p.getValue()}")

# 4. Modify a parameter (example: increase a rate constant by 2x)
target_param = model.getParameter("k3")   # parameter name varies by model
if target_param:
    old_val = target_param.getValue()
    target_param.setValue(old_val * 2.0)
    print(f"\nModified k3: {old_val} → {target_param.getValue()}")

# 5. Save modified model
writer = libsbml.SBMLWriter()
writer.writeSBMLToFile(doc, "BIOMD0000000012_modified.xml")
print("Saved modified model.")

# 6. Simulate with Tellurium (optional)
try:
    import tellurium as te
    r = te.loadSBMLModel(open("BIOMD0000000012_modified.xml").read())
    result = r.simulate(0, 100, 500)
    print(f"Simulation complete: {result.shape[0]} time points, {result.shape[1]-1} species")
    r.plot(result, title="BIOMD0000000012 modified simulation")
except ImportError:
    print("tellurium not installed — simulation step skipped")

Workflow 2: Build a Michaelis-Menten ODE Model from Scratch

Goal: Construct a full Michaelis-Menten enzyme kinetics SBML model and verify it passes validation.

import libsbml

def build_mm_model() -> libsbml.SBMLDocument:
    """Create a Michaelis-Menten enzyme kinetics SBML L3V2 model."""
    doc = libsbml.SBMLDocument(3, 2)
    model = doc.createModel()
    model.setId("michaelis_menten")
    model.setName("Michaelis-Menten Enzyme Kinetics")
    model.setTimeUnits("second")
    model.setSubstanceUnits("mole")
    model.setVolumeUnits("litre")
    model.setExtentUnits("mole")

    # Compartment
    c = model.createCompartment()
    c.setId("cell"); c.setConstant(True); c.setSize(1e-15); c.setSpatialDimensions(3)

    # Species: S (substrate), E (enzyme), ES (complex), P (product)
    for sp_id, init_conc, boundary in [
        ("S",  1e-3, False), ("E",  1e-6, False),
        ("ES", 0.0,  False), ("P",  0.0,  False)
    ]:
        sp = model.createSpecies()
        sp.setId(sp_id); sp.setCompartment("cell")
        sp.setInitialConcentration(init_conc)
        sp.setBoundaryCondition(boundary); sp.setConstant(False)
        sp.setHasOnlySubstanceUnits(False)

    # Parameters
    for p_id, p_val in [("kf", 1e6), ("kr", 1e-3), ("kcat", 0.1)]:
        p = model.createParameter()
        p.setId(p_id); p.setValue(p_val); p.setConstant(True)

    # Reactions
    def make_rxn(m, rxn_id, reacts, prods, formula):
        rxn = m.createReaction(); rxn.setId(rxn_id); rxn.setReversible(False)
        for sp, s in reacts:
            sr = rxn.createReactant(); sr.setSpecies(sp); sr.setStoichiometry(s); sr.setConstant(True)
        for sp, s in prods:
            sr = rxn.createProduct();  sr.setSpecies(sp); sr.setStoichiometry(s); sr.setConstant(True)
        kl = rxn.createKineticLaw()
        ast = libsbml.parseL3Formula(formula)
        if ast is None: raise ValueError(f"Bad formula: {formula}")
        kl.setMath(ast)

    make_rxn(model, "v_forward",  [("S",1),("E",1)], [("ES",1)], "kf * S * E * cell")
    make_rxn(model, "v_reverse",  [("ES",1)], [("S",1),("E",1)], "kr * ES * cell")
    make_rxn(model, "v_catalysis",[("ES",1)], [("P",1),("E",1)], "kcat * ES * cell")
    return doc

doc = build_mm_model()

# Validate
doc.setConsistencyChecks(libsbml.LIBSBML_CAT_UNITS_CONSISTENCY, False)  # skip units for brevity
n_errors = doc.checkConsistency()
print(f"Validation: {n_errors} issue(s)")
for i in range(n_errors):
    e = doc.getError(i)
    print(f"  [{e.getSeverityAsString()}] {e.getMessage()}")

writer = libsbml.SBMLWriter()
writer.writeSBMLToFile(doc, "michaelis_menten.xml")
print("Saved michaelis_menten.xml")
print(f"Reactions: {doc.getModel().getNumReactions()}, Species: {doc.getModel().getNumSpecies()}")

Workflow 3: Extract Stoichiometry Matrix as NumPy Array

Goal: Parse a loaded SBML model and extract the stoichiometry matrix and reaction/species lists for custom linear algebra or FBA analysis.

import libsbml
import numpy as np
import pandas as pd

doc = libsbml.readSBMLFromFile("model.xml")
model = doc.getModel()

n_species   = model.getNumSpecies()
n_reactions = model.getNumReactions()

species_ids  = [model.getSpecies(i).getId()  for i in range(n_species)]
rxn_ids      = [model.getReaction(i).getId() for i in range(n_reactions)]

# Build stoichiometry matrix S: rows=species, cols=reactions
S = np.zeros((n_species, n_reactions), dtype=float)
sp_index = {sp_id: idx for idx, sp_id in enumerate(species_ids)}

for j, rxn_id in enumerate(rxn_ids):
    rxn = model.getReaction(rxn_id)
    # Reactants: negative stoichiometry
    for k in range(rxn.getNumReactants()):
        sr = rxn.getReactant(k)
        sp_id = sr.getSpecies()
        if sp_id in sp_index:
            S[sp_index[sp_id], j] -= sr.getStoichiometry()
    # Products: positive stoichiometry
    for k in range(rxn.getNumProducts()):
        sr = rxn.getProduct(k)
        sp_id = sr.getSpecies()
        if sp_id in sp_index:
            S[sp_index[sp_id], j] += sr.getStoichiometry()

# Create a labeled DataFrame for inspection
S_df = pd.DataFrame(S, index=species_ids, columns=rxn_ids)
print("Stoichiometry matrix (S):")
print(S_df.to_string())
print(f"\nMatrix shape: {S.shape} (species × reactions)")