Document sbml functionality and add convenience functions

README.md

@@ -49,3 +49,16 @@ def get_steady_state_from_params(parameters: PyTree):
 jacobian = jax.jacrev(get_steady_state_from_params)(model.parameters)
 
 ```
+### Load a kinetic model from an sbml file
+
+```python
+from enzax.sbml import load_libsbml_model_from_url, sbml_to_enzax
+
+url = "https://raw.githubusercontent.com/dtu-qmcm/enzax/refs/heads/main/tests/data/exampleode.xml"
+
+libsbml_model = load_libsbml_model_from_url(url)
+# or to load an sbml file from your computer:
+# libsbml_model = load_libsbml_model_from_file(path_to_file)
+
+model = sbml_to_enzax(libsbml_model)
+```

docs/api/sbml.md

@@ -0,0 +1,10 @@
+#  ::: enzax.sbml
+    options:
+      show_root_heading: true
+      filters:
+        - "!check"
+      members:
+        - load_libsbml_model_from_url
+        - load_libsbml_model_from_file
+        - sbml_to_enzax
+

docs/getting_started.md

@@ -128,6 +128,30 @@ dcdt = model.dcdt(conc)
 dcdt
 ```
 
+## Load a kinetic model from SBML
+
+Enzax supports loading kinetic models from SBML files, either locally:
+
+```python
+from pathlib import Path
+from enzax.sbml import load_libsbml_model_from_file, sbml_to_enzax
+
+path = Path("path") / "to" / "sbml_file.xml"
+libsbml_model = load_libsbml_model_from_file(path)
+model = sbml_to_enzax(libsbml_model)
+```
+
+
+or from a url:
+
+```python
+from enzax.sbml import load_libsbml_model_from_url, sbml_to_enzax
+
+url = "https://raw.githubusercontent.com/dtu-qmcm/enzax/refs/heads/main/tests/data/exampleode.xml"
+libsbml_model = load_libsbml_model_from_url(url)
+model = sbml_to_enzax(libsbml_model)
+```
+
 ## Find a kinetic model's steady state
 
 Enzax provides a few example kinetic models, including [`methionine`](https://github.com/dtu-qmcm/enzax/blob/main/src/enzax/examples/methionine.py), a model of the mammalian methionine cycle.

pyproject.toml

@@ -15,6 +15,7 @@ dependencies = [
     "python-libsbml>=5.20.4",
     "sympy2jax>=0.0.5",
     "sbmlmath>=0.2.0",
+    "requests>=2.32.3",
 ]
 requires-python = ">=3.12"
 readme = "README.md"

scripts/sbml_demo.py

@@ -1,79 +1,19 @@
 import jax
 import jax.numpy as jnp
-from enzax import sbml
-from enzax.kinetic_model import KineticModelStructure, KineticModelSbml
+from enzax.sbml import load_libsbml_model, sbml_to_enzax
 from enzax.steady_state import get_kinetic_model_steady_state
 
 jax.config.update("jax_enable_x64", True)
 
-file_path = "tests/data/exampleode_names.xml"
-model_sbml = sbml.load_sbml(file_path)
-reactions_sympy = sbml.sbml_to_sympy(model_sbml)
-sym_module = sbml.sympy_to_enzax(reactions_sympy)
-
-species = [s.getId() for s in model_sbml.getListOfSpecies()]
-
-balanced_species = [
-    b.getId() for b in model_sbml.getListOfSpecies() if not b.boundary_condition
-]
-
-reactions = [reaction.getId() for reaction in model_sbml.getListOfReactions()]
-
-stoichiometry = {
-    reaction.getId(): {
-        r.getSpecies(): -r.getStoichiometry(),
-        p.getSpecies(): p.getStoichiometry(),
-    }
-    for reaction in model_sbml.getListOfReactions()
-    for r in reaction.getListOfReactants()
-    for p in reaction.getListOfProducts()
-}
-
-structure = KineticModelStructure(
-    stoichiometry=stoichiometry,
-    species=species,
-    reactions=reactions,
-    balanced_species=balanced_species,
-)
-
-parameters_local = {
-    p.getId(): p.getValue()
-    for r in model_sbml.getListOfReactions()
-    for p in r.getKineticLaw().getListOfParameters()
-}
-
-parameters_global = {
-    p.getId(): p.getValue()
-    for p in model_sbml.getListOfParameters()
-    if p.constant
-}
-
-compartments = {c.getId(): c.volume for c in model_sbml.getListOfCompartments()}
-
-unbalanced_species = {
-    u.getId(): u.getInitialConcentration()
-    for u in model_sbml.getListOfSpecies()
-    if u.boundary_condition
-}
-
-para = {
-    **parameters_local,
-    **parameters_global,
-    **compartments,
-    **unbalanced_species,
-}
-
-kinmodel_sbml = KineticModelSbml(
-    parameters=para,
-    structure=structure,
-    sym_module=sym_module,
-)
+file_path = "tests/data/brusselator.xml"
+model_libsbml = load_libsbml_model(file_path)
+model = sbml_to_enzax(model_libsbml)
 
 y0 = jnp.array([2.0, 4])
-kinmodel_sbml.flux(y0)
+model.flux(y0)
 
-kinmodel_sbml.dcdt(t=1, conc=y0)
+model.dcdt(t=1, conc=y0)
 
 guess = jnp.full((2), 0.01)
-steady_state = get_kinetic_model_steady_state(kinmodel_sbml, guess)
+steady_state = get_kinetic_model_steady_state(model, guess)
 print(steady_state)

src/enzax/sbml.py

@@ -1,11 +1,14 @@
+from pathlib import Path
+
 import libsbml
-from sbmlmath import SBMLMathMLParser
+import requests
 import sympy2jax
+from sbmlmath import SBMLMathMLParser
 
+from enzax.kinetic_model import KineticModelSbml, KineticModelStructure
 
-def load_sbml(file_path):
-    reader = libsbml.SBMLReader()
-    doc = reader.readSBML(file_path)
+
+def _get_libsbml_model_from_doc(doc):
     if doc.getModel() is None:
         raise ValueError("Failed to load the SBML model")
     elif doc.getModel().getNumFunctionDefinitions():
@@ -17,6 +20,21 @@ def load_sbml(file_path):
     return model
 
 
+def load_libsbml_model_from_file(file_path: Path) -> libsbml.Model:
+    """Load a libsbml.Model object from local file at file_path."""
+    reader = libsbml.SBMLReader()
+    doc = reader.readSBML(file_path)
+    return _get_libsbml_model_from_doc(doc)
+
+
+def load_libsbml_model_from_url(url: str) -> libsbml.Model:
+    """Load a libsbml.Model object from a url."""
+    reader = libsbml.SBMLReader()
+    with requests.get(url) as response:
+        doc = reader.readSBMLFromString(response.text)
+    return _get_libsbml_model_from_doc(doc)
+
+
 def sbml_to_sympy(model):
     reactions_sbml = model.getListOfReactions()
     reactions_sympy = [
@@ -37,3 +55,79 @@ def sbml_to_sympy(model):
 def sympy_to_enzax(reactions_sympy):
     sym_module = sympy2jax.SymbolicModule(reactions_sympy)
     return sym_module
+
+
+def get_sbml_parameters(model: libsbml.Model) -> dict:
+    kinetic_law_parameters = {
+        p.getId(): p.getValue()
+        for r in model.getListOfReactions()
+        for p in r.getKineticLaw().getListOfParameters()
+    }
+    compartment_volumes = {
+        c.getId(): c.volume for c in model.getListOfCompartments()
+    }
+    unbalanced_species = {
+        u.getId(): u.getInitialConcentration()
+        for u in model.getListOfSpecies()
+        if u.boundary_condition
+    }
+    other_parameters = {
+        p.getId(): p.getValue()
+        for p in model.getListOfParameters()
+        if p.constant
+    }
+    return {
+        **kinetic_law_parameters,
+        **compartment_volumes,
+        **unbalanced_species,
+        **other_parameters,
+    }
+
+
+def get_reaction_stoichiometry(reaction: libsbml.Reaction) -> dict[str, float]:
+    reactants = reaction.getListOfReactants()
+    products = reaction.getListOfProducts()
+    reactant_stoichiometries, product_stoichiometries = (
+        {s.getSpecies(): coeff * s.getStoichiometry() for s in list_of_species}
+        for list_of_species, coeff in [(reactants, -1.0), (products, 1.0)]
+    )
+    return reactant_stoichiometries | product_stoichiometries
+
+
+def get_sbml_structure(model_sbml: libsbml.Model) -> KineticModelStructure:
+    species = [s.getId() for s in model_sbml.getListOfSpecies()]
+    balanced_species = [
+        b.getId()
+        for b in model_sbml.getListOfSpecies()
+        if not b.boundary_condition
+    ]
+    reactions = [
+        reaction.getId() for reaction in model_sbml.getListOfReactions()
+    ]
+    stoichiometry = {
+        reaction.getId(): get_reaction_stoichiometry(reaction)
+        for reaction in model_sbml.getListOfReactions()
+    }
+    return KineticModelStructure(
+        stoichiometry=stoichiometry,
+        species=species,
+        reactions=reactions,
+        balanced_species=balanced_species,
+    )
+
+
+def get_sbml_sym_module(model: libsbml.Model):
+    reactions_sympy = sbml_to_sympy(model)
+    return sympy_to_enzax(reactions_sympy)
+
+
+def sbml_to_enzax(model: libsbml.Model) -> KineticModelSbml:
+    """Convert a KineticModelSbml object into a libsbml.Model."""
+    parameters = get_sbml_parameters(model)
+    structure = get_sbml_structure(model)
+    sym_module = get_sbml_sym_module(model)
+    return KineticModelSbml(
+        parameters=parameters,
+        structure=structure,
+        sym_module=sym_module,
+    )

tests/test_sbml.py

@@ -1,7 +1,9 @@
+from jax import numpy as jnp
 import pytest
 import importlib.resources
 from tests import data
 from enzax import sbml
+from enzax.steady_state import get_kinetic_model_steady_state
 
 brusselator_file = importlib.resources.files(data) / "brusselator.xml"
 exampleode_file = importlib.resources.files(data) / "exampleode.xml"
@@ -14,19 +16,22 @@
         exampleode_file,
     ],
 )
-def test_load_sbml(file_path):
-    sbml.load_sbml(file_path)
+def test_load_libsbml_model(file_path):
+    sbml.load_libsbml_model(file_path)
 
 
 @pytest.mark.parametrize(
-    "model",
+    ["path", "expected", "guess"],
     [
-        sbml.load_sbml(brusselator_file),
-        sbml.load_sbml(exampleode_file),
+        (
+            exampleode_file,
+            jnp.array([0.78276539, 3.65098512]),
+            jnp.array([0.01, 0.01]),
+        ),
     ],
 )
-def test_sbml_to_sympy(model):
-    sbml.sbml_to_sympy(model)
-
-
-def test_sympy_to_enzax(): ...
+def test_sbml_to_enzax(path, expected, guess):
+    libsbml_model = sbml.load_libsbml_model(path)
+    model = sbml.sbml_to_enzax(libsbml_model)
+    steady_state = get_kinetic_model_steady_state(model, guess)
+    assert jnp.isclose(steady_state, expected).all()