libsbml-network-modeling

libsbml-network-modeling

概要

libSBML は、SBML（Systems Biology Markup Language）モデルの読み込み、書き込み、作成、検証を行うためのリファレンスライブラリです。SBML は、生化学反応ネットワークを記述するためのコミュニティ標準です。ODE モデル、シグナル伝達カスケード、ゲノムスケール代謝モデルのすべてで使用されています。Python API（python-libsbml）は、コンパートメント、種、反応、速度論的法則、ルール、制約、すべての SBML 拡張を網羅する完全なオブジェクトモデルを提供します。SBML .xml ファイルとして保存されたモデルは、COPASI、Tellurium、RoadRunner、COBRApy、BioModels Database と相互運用できます。

使用する場面

• COPASI または Tellurium でシミュレーションするため、SBML 形式で新しい ODE ベースの生化学モデル（酵素反応速度論、シグナル伝達経路）を一から構築する場合
• 既存の BioModels Database モデルを読み込み、速度論的パラメーターの変更、種の追加、反応の化学量論の修正などをプログラムで行う場合
• BioModels への登録や共同研究者との共有に先立ち、SBML ファイルが仕様に準拠しているか検証する場合
• 古いシミュレーションツールとの互換性のために、SBML モデルを Level 1／2／3 間で変換する場合
• COBRApy で使用するため、FBC（Flux Balance Constraints）拡張を介してフラックスの上下限と目的関数を備えたゲノムスケール代謝モデルを構築する場合
• SBML モデルを解析し、カスタム分析用に化学量論行列、種の一覧、反応ネットワークを NumPy／pandas のデータ構造として抽出する場合
• 既に構築された代謝モデルに対して FBA、FVA、遺伝子ノックアウトを実行する必要がある場合は、代わりに cobrapy-metabolic-modeling を使用してください。libSBML は SBML ファイル自体の構築と編集に使用します
• 低レベルの XML 操作をせず、SBML の作成（Antimony 構文）と ODE シミュレーションの両方に対応する統合 Python 環境が必要な場合は、tellurium を直接使用してください

前提条件

• Python パッケージ：python-libsbml、numpy、pandas（任意。行列の抽出に使用）
• オプションのパッケージ：cobra（COBRApy。SBML の読み込み後に FBA を実行する場合）、tellurium（SBML↔Antimony の変換とシミュレーションに使用）
• データ要件：SBML ファイル（.xml）、または Python で一から構築したモデル。BioModels Database の SBML ファイルは https://www.ebi.ac.uk/biomodels/ から無料で入手できます

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

クイックスタート

SBML ファイルを読み込み、内容を調べ、パラメーター値を変更します。

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

コア API

モジュール 1：SBML の読み込みと検証

ディスクまたは文字列から SBML ファイルを読み込み、解析エラーを確認し、SBML 仕様に対する完全な検証を実行します。

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

モジュール 2：一からモデルを作成する

コンパートメント、種、反応をプログラムで追加して、完全な SBML ドキュメントを構築します。

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

モジュール 3：反応と速度論的法則の編集

化学量論と MathML の速度論的法則の式を指定して、反応を追加します。

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

モジュール 4：種とコンパートメント

種のプロパティ（初期量と濃度、境界条件、コンパートメントの体積）を確認し、変更します。

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

モジュール 5：ルールと制約

代数ルール、速度ルール、代数ルールを追加して、派生量や保存関係をモデル化します。

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()}")

モジュール 6：FBC 拡張（Flux Balance Constraints）

SBML FBC パッケージを使用し、フラックスの上下限と目的関数を備えたゲノムスケール代謝モデルを記述します。これらのモデルは COBRApy やその他の FBA ソルバーで使用できます。

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

モジュール 7：エクスポートとレベル変換

モデルをファイルまたは文字列に書き出し、SBML のレベル間で変換し、Tellurium を介して Antimony 表記にエクスポートします。

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

重要な概念

SBMLDocument、Model、プラグインアーキテクチャ

libSBML の各セッションは、Model をちょうど 1 つ所有する SBMLDocument から始まります。拡張パッケージ（FBC、qual、layout、groups、distrib）には、object.getPlugin("fbc") で取得するプラグインを介してアクセスします。プラグインを使用するには、まず doc.enablePackage(...) でドキュメント上のパッケージを有効にする必要があります。パッケージが有効になっていないドキュメントに対して getPlugin を呼び出すと、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 数式と AST

SBML では、反応速度式は MathML として保存されます。libSBML は libsbml.parseL3Formula(string) を使用して数式文字列を抽象構文木（AST）に解析し、libsbml.formulaToL3String(ast) を使用して AST を文字列に変換します。反応速度式に割り当てる前に、parseL3Formula が None 以外を返すことを必ず確認してください。None が返された場合、解析はエラーを通知せずに失敗しています。

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

よくあるワークフロー

ワークフロー 1：BioModels モデルを読み込み、パラメーターを変更してシミュレーションする

目標：BioModels モデルをダウンロードし、反応速度パラメーターを調整して、Tellurium/RoadRunner で ODE シミュレーションを実行します。

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

ワークフロー 2：Michaelis-Menten の ODE モデルを一から構築する

目標：Michaelis-Menten 酵素反応速度論の完全な SBML モデルを構築し、検証に合格することを確認します。

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()}")

ワークフロー 3：化学量論行列を NumPy 配列として抽出する

目標：読み込んだ SBML モデルを解析し、カスタム線形代数や FBA 解析用に化学量論行列と反応／種のリストを抽出します。

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

# Null-space rank as a basic model check
rank = np.linalg.matrix_rank(S)
print(f"Rank of S: {rank}")
print(f"Degrees of freedom (flux modes): {n_reactions - rank}")

ワークフロー 4：SBML FBA モデルを読み込み、COBRApy に渡す

目標：SBML FBC 形式のゲノムスケール代謝モデルを読み込み、FBA 解析のために COBRApy にロードします。

import libsbml
import cobra
import cobra.io

# Method A: use COBRApy's built-in SBML reader (wraps libSBML)
model_cobra = cobra.io.read_sbml_model("iJO1366.xml")
print(f"COBRApy model: {model_cobra.id}")
print(f"  Reactions: {len(model_cobra.reactions)}")
print(f"  Metabolites: {len(model_cobra.metabolites)}")
print(f"  Genes: {len(model_cobra.genes)}")

# Run FBA
solution = model_cobra.optimize()
print(f"\nFBA objective value: {solution.objective_value:.4f}")
print(f"Status: {solution.status}")

# Method B: use libSBML to inspect FBC metadata before loading into COBRApy
doc = libsbml.readSBMLFromFile("iJO1366.xml")
model = doc.getModel()
fbc = model.getPlugin("fbc")

if fbc:
    n_obj = fbc.getNumObjectives()
    active_obj_id = fbc.getActiveObjectiveId()
    print(f"\nlibSBML FBC: {n_obj} objective(s), active='{active_obj_id}'")
    obj = fbc.getObjective(active_obj_id)
    if obj:
        print(f"Objective type: {obj.getType()}")
        for i in range(obj.getNumFluxObjectives()):
            fo = obj.getFluxObjective(i)
            print(f"  {fo.getReaction()} (coeff={fo.getCoefficient()})")

主なパラメーター

パラメーター	モジュール	デフォルト	範囲／オプション	効果
level	SBMLDocument	3	1、2、3	SBML レベルです。新しいモデルでは L3 を、従来のツールとの互換性には L2 を使用します
version	SBMLDocument	2	L3 では 1～2、L2 では 1～4	レベル内の SBML バージョンです。L3V2 が現在の標準です
initialConcentration	Species	0.0	任意の浮動小数点数	初期モル濃度です。initialAmount とは相互に排他的です
hasOnlySubstanceUnits	Species	False	True、False	True の場合、反応速度式は物質量（mol）を参照し、False の場合は濃度（M）を参照します
boundaryCondition	Species	False	True、False	True の場合、ODE ソルバーはこの種の量を変更しません。外部入力に使用します
constant	Species/Parameter	True (Param)	True、False	代入規則の対象および変更可能なパラメーターには False が必要です
strict	FBC plugin	True	True、False	FBC 厳密モードでは、すべてのフラックス境界値をパラメーターとして定義する必要があります
targetLevel / targetVersion	ConversionProperties	—	L/V の整数	doc.convert() によるレベル／バージョン変換のターゲットです
LIBSBML_CAT_UNITS_CONSISTENCY	checkConsistency	True	True、False	検証時の単位次元チェックを有効／無効にします

ベストプラクティス

1. parseL3Formula の戻り値は必ず確認してください：この関数は、不正な入力に対して例外を発生させずに None を返します。None を kl.setMath() に割り当てると、構文解析は通るものの、検証に失敗する反応速度式が欠落したモデルが作成されます。

   ast = libsbml.parseL3Formula("Vmax * S / (Km + S)")
   if ast is None:
       raise ValueError("Formula parse failed")
   kl.setMath(ast)
   

2. 代入ルールの対象には constant=False を設定してください：AssignmentRule または RateRule の対象となるパラメーターや種には、constant を False に設定する必要があります。True にすると制約違反が生じ、整合性チェックに失敗します。

3. 反応速度式では反応速度に区画の体積を掛けてください：SBML の反応進行量の単位はモル（または分子数）であるため、反応速度の単位は反応進行量／時間でなければなりません。濃度で表される種については、区画の大きさを掛けます：kf * S * E * compartment_volume。この係数を省略することが、反応速度式で最もよくある単位の誤りです。

4. 使用するパッケージだけを有効にしてください：不要な拡張機能（layout、groups）に対して doc.enablePackage() を呼び出すと、名前空間宣言が追加され、一部の後続ツールで問題が発生します。FBC は FBA／FVA モデルでのみ有効にし、純粋な ODE モデルでは無効のままにしてください。

5. 手早く読み込むには readSBMLFromFile（トップレベル関数）を使用してください：便利な関数 libsbml.readSBMLFromFile(path) は、SBMLReader インスタンスを作成してそのインスタンスの readSBMLFromFile を呼び出す方法と同等です。どちらも SBMLDocument を返すため、記述が簡潔な方を選んでください。

6. 保存前と変換後に検証してください：writeSBMLToFile を呼び出す直前に doc.checkConsistency() を実行し、さらにレベル／バージョン変換の後にも実行してください。特に単位に関して、変換によって新たな警告が発生することがあります。

よくあるレシピ

レシピ：すべての反応と反応速度式を一覧表示する

使用する場面：SBML モデルを変更する前に監査し、すべての反応とその速度則を記録します。

import libsbml

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

print(f"{'Reaction':<20} {'Reversible':<12} {'Kinetic Law'}")
print("-" * 80)
for i in range(model.getNumReactions()):
    rxn = model.getReaction(i)
    if rxn.isSetKineticLaw():
        formula = libsbml.formulaToL3String(rxn.getKineticLaw().getMath())
    else:
        formula = "(no kinetic law)"
    rev = "reversible" if rxn.getReversible() else "irreversible"
    print(f"{rxn.getId():<20} {rev:<12} {formula}")

レシピ：複数のパラメーターを一括更新する

使用する場面：感度分析を行うために、一連の速度定数をさまざまな値に変更し、値ごとに SBML モデルを再保存します。

import libsbml
import copy

doc = libsbml.readSBMLFromFile("model.xml")
writer = libsbml.SBMLWriter()

# Parameter sweep: vary kcat and Km
sweep = [
    {"kcat": 0.05, "Km": 0.01},
    {"kcat": 0.10, "Km": 0.01},
    {"kcat": 0.20, "Km": 0.01},
    {"kcat": 0.10, "Km": 0.05},
]

for idx, params in enumerate(sweep):
    # Re-read fresh copy each iteration to avoid cumulative edits
    doc_i = libsbml.readSBMLFromFile("model.xml")
    model_i = doc_i.getModel()
    for p_id, p_val in params.items():
        p = model_i.getParameter(p_id)
        if p:
            p.setValue(p_val)
    fname = f"model_sweep_{idx:03d}.xml"
    writer.writeSBMLToFile(doc_i, fname)
    print(f"Saved {fname}: {params}")

レシピ：すべての種の初期条件を Dict として抽出する

使用する場面：SBML で定義された初期条件を使って、カスタム ODE ソルバー（例：scipy.integrate.solve_ivp）を初期化します。

import libsbml

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

initial_conditions = {}
for i in range(model.getNumSpecies()):
    sp = model.getSpecies(i)
    if sp.isSetInitialConcentration():
        initial_conditions[sp.getId()] = sp.getInitialConcentration()
    elif sp.isSetInitialAmount():
        # Convert amount to concentration using compartment volume
        comp = model.getCompartment(sp.getCompartment())
        vol = comp.getSize() if comp and comp.isSetSize() else 1.0
        initial_conditions[sp.getId()] = sp.getInitialAmount() / vol
    else:
        initial_conditions[sp.getId()] = 0.0

print("Initial conditions (concentration in model units):")
for sp_id, val in initial_conditions.items():
    print(f"  {sp_id}: {val:.6g}")

トラブルシューティング

問題	原因	解決策
ImportError: No module named 'libsbml'	パッケージがインストールされていません	pip install python-libsbml を実行します。インポート名は python_libsbml ではなく libsbml であることに注意してください
parseL3Formula が None を返します	数式文字列の形式が不正です（演算子が誤っている、関数が未定義など）	数式の構文を確認します。正常な AST に対して libsbml.formulaToL3String を使うと、期待される形式を確認できます。* は乗算、^ または pow() は累乗です
整合性チェックで単位エラーが報告されます	反応速度式にコンパートメント体積の係数がありません	反応速度式にコンパートメント体積を乗算します：kf * S * E * V。モデルで substance_units = "mole" と volume_units = "litre" を設定します
getPlugin("fbc") が None を返します	ドキュメントで FBC パッケージが有効になっていません	モデルの読み取りまたは構築の前に doc.enablePackage(libsbml.FbcExtension.getXmlnsL3V1V2(), "fbc", True) を呼び出します
レベル／バージョン変換がゼロ以外のコードを返します	ソースモデルに、変換先のレベルでサポートされていない機能が含まれています	変換後に doc.getNumErrors() を確認します。SBML L1 には重大な制限があります（コンパートメントや単位がありません）。最低限の変換先として L2V4 を推奨します
割り当てルールの対象が「model is overdetermined」エラーを引き起こします	種またはパラメーターが constant=True であるにもかかわらず、ルールの対象になっています	ルールの対象で constant=False を設定します。constant=True は値が固定され、ルールで上書きできないことを意味します
カスタム構築した SBML で COBRApy の read_sbml_model が失敗します	FBC の strict=True が設定されていますが、フラックス境界がパラメーターとして定義されていません	すべての反応の FBC プラグインで、setLowerFluxBound と setUpperFluxBound が既存のパラメーター ID を参照していることを確認します
大規模モデルの読み込みに時間がかかります（10 秒超）	SBML ファイルが非常に大きいです（ゲノムスケールモデル、反応数 1 万以上）	正常な動作です。libSBML の XML 解析はシングルスレッドです。readSBMLFromFile を使用し（文字列ベースの読み込みは避け）、ループ内での再読み込みを避けます

関連スキル

• cobrapy-metabolic-modeling — SBML／JSON 代謝モデルに対する FBA、FVA、遺伝子ノックアウト、フラックスサンプリング。libSBML でモデルを構築または編集し、COBRApy で解析します
• string-database-ppi — タンパク質間相互作用ネットワーク。PPI データを SBML qual 拡張にエクスポートし、論理ネットワークモデルを作成します
• reactome-database — 経路データソース。Reactome はヒト経路の SBML エクスポートを提供しており、libSBML で読み込んで解析できます
• brenda-database — libSBML の反応速度式に実験測定値を設定するための速度論的パラメーター（Km、Vmax、kcat）のデータソース
• networkx-graph-analysis — SBML から反応ネットワークを抽出した後に、グラフのトポロジー（最短経路、連結性、中心性）を解析するために使用します
• sympy-symbolic-math — SBML から解析した反応速度式の記号操作。formulaToL3String と組み合わせて、解析的な定常状態の導出に使用します

参考資料

• libSBML Python ドキュメント — Python API の完全なリファレンス
• SBML 仕様ポータル — SBML Level 3 コアおよびパッケージ仕様の公式資料
• libSBML GitHub リポジトリ — ソースコード、課題トラッカー、リリースノート
• BioModels Database — キュレーション済みの SBML モデルリポジトリ（手作業でキュレーションされたモデルが 1,000 件以上）
• Hucka et al. (2003) Bioinformatics 19(4):524–531 — SBML 仕様の原著論文
• Keating et al. (2020) Molecular Systems Biology 16(8):e9110 — SBML Level 3 と現在の拡張エコシステムの概要