added adaptor method and one test

src/emll/test_antimony_utils.py

@@ -0,0 +1,41 @@
+"""Test file for Antimony to Cobra conversion"""
+
+import logging
+from pathlib import Path
+
+import cobra
+import pytest
+import tellurium as te
+
+from emll.util import ant_to_cobra
+
+logging.getLogger("cobra").setLevel(logging.ERROR)
+
+HERE = Path(__file__).parent.resolve()
+MODEL = HERE.joinpath("test_models/unlabeled18.ant")
+
+
+@pytest.mark.parametrize(
+    "antimony_path",
+    [
+        str(MODEL),
+    ],
+)
+def test_antimony_to_cobra_conversion(antimony_path):
+    """
+    Compares the stoichiometric matrices of original antimony file
+    and cobra-compatible sbml file to see if they are the same.
+
+    input: path to original antimony file (non-cobra-compatible)
+    """
+    ant_to_cobra(antimony_path)
+
+    r = te.loada(antimony_path)
+    N_te = r.getFullStoichiometryMatrix()
+
+    file_name = antimony_path.split(".")[0]
+    model = cobra.io.read_sbml_model(file_name + "_cobra.xml")
+    N_cobra = cobra.util.create_stoichiometric_matrix(model)
+
+    assert N_te.shape == N_cobra.shape
+    assert (N_te == N_cobra).all()

src/emll/test_models/unlabeled18.ant

@@ -0,0 +1,148 @@
+
+v2: $A -> B; e_v2*(Vm2 / Km2 *(A-B/Keq2)) / (1 + A/Km2 + B/Km3)
+v3: B -> R; e_v3*(Vm3 / Km4 *(B-R/Keq3)) / (1 + B/Km4 + R/Km5)
+v13: R -> ; e_v13*(Vm13 / Km22 *(R/Keq13)) / (1 + R/Km22)
+
+v4: B -> C; e_v4*(Vm4 / Km6 *(B-C/Keq4)) / (1 + B/Km6 + C/Km7)
+v5: C -> D; e_v5*(Vm5 / Km7 *(C-D/Keq5)) / (1 + C/Km7 + D/Km8)
+v6: D -> E; e_v6*(Vm6 / Km9 *(D-E/Keq6)) / (1 + D/Km9 + E/Km10)
+v7: E ->; e_v7*(Vm7 / Km11 *(E/Keq7)) / (1 + E/Km11)
+v8: D -> F; e_v8*(Vm8 / Km12 *(D-F/Keq8)) / (1 + D/Km12 + F/Km13)
+v9: D -> G; e_v9*(Vm9 / Km14 *(D-G/Keq9)) / (1 + D/Km14 + G/Km15)
+v10: G -> H; e_v10*(Vm10 / Km16 *(G-H/Keq10)) / (1 + G/Km16 + H/Km17)
+v11: H -> I; e_v11*(Vm11 / Km18 *(H-I/Keq11)) / (1 + H/Km18 + I/Km19)
+v12: I -> $J; e_v12*(Vm12 / Km20 *(I-J/Keq12)) / (1 + I/Km20 + J/Km21)
+
+v14: F -> ; e_v14*(Vm14 / Km23 *(F/Keq14)) / (1 + F/Km23)
+v15: F -> L; e_v15*(Vm15 / Km24 *(F-L/Keq15)) / (1 + F/Km24 + L/Km25)
+v16: L -> M; e_v16*(Vm16 / Km26 *(L-M/Keq16)) / (1 + L/Km26 + M/Km27)
+v17: M -> N; e_v17*(Vm17 / Km28 *(M-N/Keq17)) / (1 + M/Km28 + N/Km29)
+v18: N-> O; e_v18*(Vm18 / Km30 *(N-O/Keq18)) / (1 + N/Km30 + O/Km31)
+v19: O-> $P; e_v19*(Vm19 / Km31 *(O-P/Keq19)) / (1 + O/Km31 + P/Km32)
+
+v21: O -> $Q; e_v21*(Vm21 / Km34 *(O-Q/Keq21)) / (1 + O/Km34 + Q/Km35)
+
+e_v1 = 1;
+e_v2 = 1;
+e_v3 = 1;
+e_v4 = 1;
+e_v5 = 1;
+e_v6 = 1;
+e_v7 = 1;
+e_v8 = 1;
+e_v9 = 1;
+e_v10 = 1;
+e_v11 = 1;
+e_v12 = 1;
+e_v13 = 1;
+e_v14 = 1;
+e_v15 = 1;
+e_v16 = 1;
+e_v17 = 1;
+e_v18 = 1;
+e_v19 = 1;
+e_v20 = 1;
+e_v21 = 1;
+e_v22 = 1;
+
+Vm1 = 1.448;
+Vm2 = 4.881;
+Vm3 = 3.556;
+Vm4 = 9.404;
+Vm5 = 7.653;
+Vm6 = 7.487;
+Vm7 = 9.037;
+Vm8 = 0.834;
+Vm9 = 5.522;
+Vm10 = 5.845;
+Vm11 = 9.619;
+Vm12 = 2.921;
+Vm13 = 2.408;
+Vm14 = 1.003;
+Vm15 = 0.164;
+Vm16 = 9.295;
+Vm17 = 6.699;
+Vm18 = 7.852;
+Vm19 = 2.817;
+Vm20 = 5.864;
+Vm21 = 0.64;
+Vm22 = 4.856;
+
+Keq1 = 0.6875;
+Keq2 = 2.155;
+Keq3 = 9.474;
+Keq4 = 7.309;
+Keq5 = 2.539;
+Keq6 = 2.133;
+Keq7 = 5.182;
+Keq8 = 0.257;
+Keq9 = 2.075;
+Keq10 = 4.247;
+Keq11 = 3.742;
+Keq12 = 4.636;
+Keq13 = 2.776;
+Keq14 = 5.868;
+Keq15 = 8.639;
+Keq16 = 1.175;
+Keq17 = 5.174;
+Keq18 = 1.321;
+Keq19 = 7.169;
+Keq20 = 3.961;
+Keq21 = 5.654;
+Keq22 = 1.833;
+
+Km1 = 0.398;
+Km2 = 0.21;
+Km3 = 0.186;
+Km4 = 0.944;
+Km5 = 0.74;
+Km6 = 0.49;
+Km7 = 0.227;
+Km8 = 0.254;
+Km9 = 0.058;
+Km10 = 0.434;
+Km11 = 0.312;
+Km12 = 0.696;
+Km13 = 0.378;
+Km14 = 0.18;
+Km15 = 0.025;
+Km16 = 0.067;
+Km17 = 0.679;
+Km18 = 0.454;
+Km19 = 0.537;
+Km20 = 0.897;
+Km21 = 0.99;
+Km22 = 0.217;
+Km23 = 0.663;
+Km24 = 0.263;
+Km25 = 0.021;
+Km26 = 0.758;
+Km27 = 0.32;
+Km28 = 0.383;
+Km29 = 0.588;
+Km30 = 0.831;
+Km31 = 0.629;
+Km32 = 0.873;
+Km33 = 0.274;
+Km34 = 0.798;
+Km35 = 0.186;
+Km36 = 0.953;
+
+A = 10; 
+B = 0.19;
+C = 0.511;
+D = 0.224;
+E = 0.098;
+F = 0.862;
+G = 0.973;
+H = 0.961;
+I = 0.907;
+J = 3.32;
+# K = 0.333;
+L = 0.081;
+M = 0.407;
+N = 0.232;
+O = 0.132;
+P = 0.53;
+Q = 2.26;
+R = 0.011;

src/emll/util.py

@@ -2,6 +2,12 @@
 import scipy as sp
 import pytensor.tensor as at
 import pymc as pm
+from pathlib import Path
+
+import tellurium as te
+import libsbml
+
+import os
 
 
 def create_elasticity_matrix(model):
@@ -25,11 +31,15 @@ def create_elasticity_matrix(model):
                 array[r_ind(reaction), m_ind(metabolite)] = -np.sign(stoich)
 
             # Irrevesible in forward direction, only assign if met is reactant
-            elif (not reaction.reversibility) & (reaction.upper_bound > 0) & (stoich < 0):
+            elif (
+                (not reaction.reversibility) & (reaction.upper_bound > 0) & (stoich < 0)
+            ):
                 array[r_ind(reaction), m_ind(metabolite)] = -np.sign(stoich)
 
             # Irreversible in reverse direction, only assign in met is product
-            elif (not reaction.reversibility) & (reaction.lower_bound < 0) & (stoich > 0):
+            elif (
+                (not reaction.reversibility) & (reaction.lower_bound < 0) & (stoich > 0)
+            ):
                 array[r_ind(reaction), m_ind(metabolite)] = -np.sign(stoich)
 
     return array
@@ -43,7 +53,8 @@ def create_Ey_matrix(model):
 
     boundary_indexes = [model.reactions.index(r) for r in model.medium.keys()]
     boundary_directions = [
-        1 if r.products else -1 for r in model.reactions.query(lambda x: x.boundary, None)
+        1 if r.products else -1
+        for r in model.reactions.query(lambda x: x.boundary, None)
     ]
     ny = len(boundary_indexes)
     Ey = np.zeros((len(model.reactions), ny))
@@ -151,7 +162,9 @@ def initialize_elasticity(
         name = "ex"
 
     if m_compartments is not None:
-        assert r_compartments is not None, "reaction and metabolite compartments must both be given"
+        assert (
+            r_compartments is not None
+        ), "reaction and metabolite compartments must both be given"
 
         regulation_array = np.array(
             [[a in b for a in m_compartments] for b in r_compartments]
@@ -220,3 +233,111 @@ def initialize_elasticity(
     E = flat_e_entries[flat_indexer].reshape(N.T.shape)
 
     return E
+
+
+def ant_to_cobra(antimony_path):
+    """
+    This method takes in an Antimony file and converts it to a Cobra-
+    friendly format by removing all boundary species and replacing all
+    rate laws with a constant. Both an Antimony file and an SBML are
+    produced, and the file name is returned.
+    """
+    output_name = antimony_path.split("/")[-1]
+    output_name = output_name.split(".")[0]
+    output_path = Path(antimony_path).parent
+
+    # load normal Antimony file
+    with open(antimony_path, "r") as file:
+        # reprint the antimony model into a temp file to ensure proper
+        # headers
+        # load antimony file
+        r = te.loada(antimony_path)
+        # print antimony file to temp
+        with open("tempA7K8L2P4W9.txt", "w") as f:
+            f.write(r.getCurrentAntimony())
+        old_ant = "tempA7K8L2P4W9.txt"
+
+    with open(old_ant, "r") as file:
+        lines = file.readlines()
+        section_indices = []
+        c_index = -1
+        for i, line in enumerate(lines):
+            if "//" in line:
+                section_indices.append(i)
+            if "// Compartment" in line and c_index == -1:
+                c_index = i
+        next_section = section_indices.index(c_index) + 1
+
+        with open(old_ant, "r") as file:
+            lines = file.readlines()[c_index : section_indices[next_section]]
+
+        with open(f"{output_path}/{output_name}_cobra.ant", "w") as f:
+            f.write("")
+
+        for line in lines:
+            line = line.strip()
+            if "$" not in line:
+                with open(f"{output_path}/{output_name}_cobra.ant", "a") as f:
+                    f.write(line + "\n")
+            else:
+                no_bd_sp = [i for i in line.split(",") if "$" not in i]
+                no_bd_sp = [i for i in no_bd_sp if i != ""]
+
+                line = ",".join(no_bd_sp)
+                if line != "" and line[-1] != ";":
+                    line += ";"
+                if "species" not in line:
+                    line = "species" + line
+                with open(f"{output_path}/{output_name}_cobra.ant", "a") as f:
+                    f.write(line + "\n")
+
+    r = te.loada(antimony_path)
+    doc = libsbml.readSBMLFromString(r.getSBML())
+    model = doc.getModel()
+
+    bd_sp = r.getBoundarySpeciesIds()
+
+    reactants_list = []
+    products_list = []
+
+    for n in range(len(r.getReactionIds())):
+        rxn = model.getReaction(n)
+        reactants = []
+        products = []
+
+        for reactant in range(rxn.getNumReactants()):
+            stoich = rxn.getReactant(reactant).getStoichiometry()
+            if stoich == 1:
+                reactants.append(rxn.getReactant(reactant).species)
+            else:
+                reactants.append(str(stoich) + " " + rxn.getReactant(reactant).species)
+        reactants_list.append([i for i in reactants if i not in bd_sp])
+
+        for product in range(rxn.getNumProducts()):
+            stoich = rxn.getProduct(product).getStoichiometry()
+            if stoich == 1:
+                products.append(rxn.getProduct(product).species)
+            else:
+                products.append(str(stoich) + " " + rxn.getProduct(product).species)
+        products_list.append([i for i in products if i not in bd_sp])
+
+    for i in range(len(reactants_list)):
+        r1 = " + ".join(reactants_list[i])
+        p1 = " + ".join(products_list[i])
+        with open(f"{output_path}/{output_name}_cobra.ant", "a") as f:
+            f.write(r.getReactionIds()[i] + ": " + r1 + " -> " + p1 + "; 1;\n")
+
+    with open(f"{output_path}/{output_name}_cobra.ant", "a") as f:
+        f.write("\n")
+
+    for sp in r.getFloatingSpeciesIds():
+        with open(f"{output_path}/{output_name}_cobra.ant", "a") as f:
+            f.write(sp + " = 1;\n")
+
+    with open(f"{output_path}/{output_name}_cobra.xml", "w") as f:
+        f.write(te.loada(f"{output_path}/{output_name}_cobra.ant").getCurrentSBML())
+
+    if "tempA7K8L2P4W9.txt" in os.listdir():
+        os.remove("tempA7K8L2P4W9.txt")
+
+    return f"{output_path}/{output_name}_cobra"