libsbml-network-modeling

libsbml-network-modeling

概述

libSBML 是用于读取、编写、创建和验证 SBML（系统生物学标记语言）模型的参考库。SBML 是编码生化反应网络的社区标准——ODE 模型、信号级联和全基因组规模代谢模型都使用它。Python API（python-libsbml）提供了完整的对象模型，涵盖区室、物种、反应、动力学定律、规则、约束以及所有 SBML 扩展。保存为 SBML .xml 文件的模型可与 COPASI、Tellurium、RoadRunner、COBRApy 和 BioModels Database 互操作。

适用场景

• 从头开始以 SBML 格式构建新的基于 ODE 的生化模型（酶动力学、信号通路），用于在 COPASI 或 Tellurium 中进行模拟
• 读取并以编程方式修改现有的 BioModels Database 模型——更改动力学参数、添加物种或修补反应化学计量
• 在提交到 BioModels 或与协作者共享之前，根据规范验证 SBML 文件
• 在 Level 1/2/3 之间转换 SBML 模型，以兼容较旧的模拟工具
• 通过 FBC（通量平衡约束）扩展构建具有通量边界和目标函数的全基因组规模代谢模型，以便与 COBRApy 配合使用
• 解析 SBML 模型，将化学计量矩阵、物种列表或反应网络提取为 NumPy/pandas 数据结构，以便进行自定义分析
• 如果需要对已构建的代谢模型运行 FBA、FVA 或基因敲除，请改用 cobrapy-metabolic-modeling——libSBML 用于构建和编辑 SBML 文件本身
• 如果需要集成式 Python 环境，以便同时使用 Antimony 语法编写 SBML 并进行 ODE 模拟，而不必操作底层 XML，请直接使用 tellurium

先决条件

• Python 包：python-libsbml、numpy、pandas（可选，用于提取矩阵）
• 可选包：cobra（COBRApy，用于加载 SBML 后运行 FBA）、tellurium（用于 SBML↔Antimony 转换和模拟）
• 数据要求：SBML 文件（.xml），或在 Python 中从头构建；可从 https://www.ebi.ac.uk/biomodels/ 免费获取 BioModels Database 的 SBML 文件

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 扩展（通量平衡约束）

使用 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 会话都从一个 SBMLDocument 开始，该文档恰好拥有一个 Model。扩展包（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 转换回字符串。将 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 格式的基因组规模代谢模型，并将其加载到 COBRApy 中进行 FBA 分析。

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 Level；所有新模型均使用 L3；为兼容旧版工具使用 L2
version	SBMLDocument	2	L3：1–2；L2：1–4	Level 内的 SBML Version；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 插件	True	True、False	FBC 严格模式要求所有通量边界均定义为参数
targetLevel / targetVersion	ConversionProperties	—	L/V 整数	doc.convert() 进行 Level/Version 转换时的目标
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. 只启用实际使用的包：对不必要的扩展（布局、分组）调用 doc.enablePackage() 会添加命名空间声明，使某些下游工具难以处理。仅对 FBA/FVA 模型启用 FBC；纯 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；注意，导入名称是 libsbml，而不是 python_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)
Level/版本转换返回非零代码	源模型包含目标级别不支持的特性	转换后检查 doc.getNumErrors()；SBML L1 有严重限制（不支持区室和单位）；建议至少使用 L2V4 作为目标
赋值规则目标导致“模型超定”错误	物种或参数为 constant=True，但被规则指定为目标	对规则目标设置 constant=False；constant=True 表示该值固定，不能被规则覆盖
COBRApy 的 read_sbml_model 无法读取自定义构建的 SBML	FBC 的 strict=True，但通量边界未定义为参数	确保每个反应的 FBC 插件都有 setLowerFluxBound 和 setUpperFluxBound，且它们指向已存在的参数 ID
读取大型模型耗时较长（>10 秒）	SBML 文件非常大（基因组规模模型，含 10k+ 个反应）	这是正常现象——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 — 动力学参数来源（Km、Vmax、kcat），用于将实验测量值填入 libSBML 动力学定律
• 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 等（2003）Bioinformatics 19(4):524–531 — SBML 原始规范论文
• Keating 等（2020）Molecular Systems Biology 16(8):e9110 — SBML Level 3 及当前扩展生态概述