Merge pull request #25 from biocircuits/Saving-SBML

Fixed more SBML loading and writing bugs
biocircuits · Apr 23, 2020 · 2261662 · 2261662
2 parents 86bccca + 93200ad
commit 2261662
Show file tree

Hide file tree

Showing 5 changed files with 341 additions and 269 deletions.
diff --git a/bioscrape/sbmlutil.py b/bioscrape/sbmlutil.py
@@ -97,23 +97,27 @@ def import_sbml(sbml_file, bioscrape_model = None, input_printout = False):
             reactantspecies = model.getSpecies(reactant.getSpecies())
             reactantspecies_id = reactantspecies.getId()
             if reactantspecies_id in allspecies:
-                reactant_list.append(reactantspecies_id)
+                for i in range(int(reactant.getStoichiometry())):
+                    reactant_list.append(reactantspecies_id)
             else:
                 warnings.warn('Reactant in reaction {0} not found in the list of species in the SBML model.'
                  + ' The species will be added with zero initial amount'.format(reaction.getId()))
                 allspecies[reactantspecies_id] = 0.0
-                reactant_list.append(reactantspecies_id)
+                for i in range(int(reactant.getStoichiometry())):
+                    reactant_list.append(reactantspecies_id)
 
         for product in reaction.getListOfProducts():
             productspecies = model.getSpecies(product.getSpecies())
             productspecies_id = productspecies.getId()
             if productspecies_id in allspecies:
-                product_list.append(productspecies_id)
+                for i in range(int(product.getStoichiometry())):
+                    product_list.append(productspecies_id)
             else:
                 warnings.warn('Reactant in reaction {0} not found in the list of species in the SBML model.' + 
                             ' The species will be added with zero initial amount'.format(reaction.getId()))
                 allspecies[productspecies_id] = 0.0
-                product_list.append(productspecies_id)
+                for i in range(int(product.getStoichiometry())):
+                    product_list.append(productspecies_id)
 
         #Identify propensities based upon annotations
         annotation_string = reaction.getAnnotationString()
@@ -240,6 +244,21 @@ def _add_underscore_to_parameters(formula, parameters):
 
     return str(sympy_rate)
 
+def _get_species_list_in_formula(formula, species):
+    sympy_rate = sympy.sympify(formula, _clash1)
+    nodes = [sympy_rate]
+    index = 0
+    while index < len(nodes):
+        node = nodes[index]
+        index += 1
+        nodes.extend(node.args)
+    species_return = []
+    for node in nodes:
+        if type(node) == sympy.Symbol:
+            if node.name in species:
+                species_return.append(node.name)
+    return species_return
+
 def _remove_underscore_from_parameters(formula, parameters):
     sympy_rate = sympy.sympify(formula, _clash1)
     nodes = [sympy_rate]
@@ -309,8 +328,6 @@ def add_species(model, compartment, species, debug=False, initial_concentration=
     species_name = species
 
     # Construct the species ID 
-    # NOTE: These changes to the species id are not propagated into KineticLaw objects etc.
-    # TODO: Fix this. 
     species_id = species_sbml_id(species_name, model.getSBMLDocument())
 
     if debug: print("Adding species", species_name, species_id)
@@ -357,23 +374,19 @@ def add_rule(model, rule_id, rule_type, rule_variable, rule_formula, **kwargs):
         # Simply create SBML assignment rule type. For additive rule type as well, 
         # AssignmentRule type of SBML will work as $s_0$ is the artificial species that 
         # exists in the bioscrape model. 
+        if rule_variable[0] == '_':
+            rule_variable = rule_variable.replace('_','',1)
         for param in model.getListOfParameters():
             if rule_variable == param.getId():
                 param.setConstant(False)
-        # allparams = {}
-        # for p in model.getListOfParameters():
-        #     pid = p.getId()
-        #     pid = '_' + pid
-        #     allparams[pid] = p.getValue()
-        # rule_formula = _remove_underscore_from_parameters(rule_formula, allparams)
-
-        # General propensity add modifier
-        # Rules params constant false
-        # Rules underscore thing
+        allparams = {}
+        for p in model.getListOfParameters():
+            pid = p.getId()
+            pid = '_' + pid
+            allparams[pid] = p.getValue()
+        rule_formula = _remove_underscore_from_parameters(rule_formula, allparams)
         rule = model.createAssignmentRule()
         rule.setId(rule_id)
-        if rule_variable[0] == '_':
-            rule_variable = rule_variable.replace('_','',1)
         rule.setVariable(rule_variable)
         rule.setFormula(rule_formula)  
     return rule 
@@ -400,7 +413,11 @@ def add_reaction(model, inputs_list, outputs_list,
     reaction.setId(trans.getValidIdForName(reaction_id))
     ratestring = "" #Stores the string representing the rate function
     annotation_dict = {"type":propensity_type}
-
+    allspecies = []
+    for s in model.getListOfSpecies():
+        sid = s.getId()
+        allspecies.append(sid)
+
     allparams = {}
     for p in model.getListOfParameters():
         pid = p.getId()
@@ -420,7 +437,6 @@ def add_reaction(model, inputs_list, outputs_list,
         annotation_dict["n"] = propensity_params['n']
 
     elif propensity_type == "general":
-        ratestring = propensity_params['rate']
         annotation_dict["rate"] = propensity_params['rate']
     else:
         raise ValueError(propensity_type+" is not a supported propensity_type")
@@ -547,7 +563,14 @@ def add_reaction(model, inputs_list, outputs_list,
 
         annotation_dict["s1"] = s_species_id
         annotation_dict["d"] = d_species_id
-
+    elif propensity_type == "general":
+        ratestring = propensity_params['rate']
+        species_list = _get_species_list_in_formula(ratestring, allspecies)
+        for s in species_list:
+            if s not in reactants_list and s not in products_list:
+                modifier = reaction.createModifier()
+                modifier.setSpecies(s)
+
     ratestring = _remove_underscore_from_parameters(ratestring, allparams)
     # Set the ratelaw to the ratestring
     math_ast = libsbml.parseL3Formula(ratestring)
@@ -783,21 +806,6 @@ def getValidIdForName(self, name):
     #  #
 
 
-#  # Returns a list of all ids from the given list of elements
-#  #
-def getAllIds(allElements):
-    result = []
-    if (allElements == None or allElements.getSize() == 0):
-        return result
-
-    for i in range(0, allElements.getSize()):
-        current = allElements.get(i)
-        if (current.isSetId() \
-           and current.getTypeCode() != libsbml.SBML_LOCAL_PARAMETER):
-            result.append(current.getId())
-    return result
-
-
 def getSpeciesByName(model, name, compartment=''):
     '''
     Returns a list of species in the Model with the given name
@@ -824,5 +832,5 @@ def getSpeciesByName(model, name, compartment=''):
     elif not species_found:
         raise ValueError('The species ' + name + ' not found.')
     else:
-        warn('Multiple species with name ' + name + ' found. Returning a list')
+        warnings.warn('Multiple species with name ' + name + ' found. Returning a list')
         return species_found
diff --git a/bioscrape/types.pyx b/bioscrape/types.pyx
@@ -2234,7 +2234,7 @@ cdef class Model:
             add_parameter(model = model, param_name=p, param_value = val)
 
         for s in self.get_species():
-            add_species(model=model, compartment=model.getCompartment(0), species=s, initial_concentration=self.get_species_value(s))
+            add_species(model = model, compartment=model.getCompartment(0), species=s, initial_concentration=self.get_species_value(s))
 
         rxn_count = 0
         for rxn_tuple in self.reaction_definitions:
@@ -2254,7 +2254,7 @@ cdef class Model:
             # Extract the rule variable id from rule_dict:
             equation = rule_dict['equation']
             split_eqn = [s.strip() for s in equation.split('=') ]
-            assert(len(split_eqn) == 2)
+            assert(len(split_eqn) == 2) # Checking rule_dict equation structure.
             # Extract the rule formula for the variable above from rule_dict:
             rule_formula = split_eqn[1]
             rule_variable = split_eqn[0]

diff --git a/examples/Basic Examples - START HERE.ipynb b/examples/Basic Examples - START HERE.ipynb
diff --git a/examples/models/represillator_sbml_back.xml b/examples/models/represillator_sbml_back.xml
@@ -20,29 +20,29 @@
       <species id="Z" name="Z" compartment="default" initialConcentration="0" hasOnlySubstanceUnits="false" boundaryCondition="false" constant="false"/>
     </listOfSpecies>
     <listOfParameters>
-      <parameter id="beta" value="0.2" constant="true"/>
-      <parameter id="alpha0" value="0.216404256133345" constant="true"/>
-      <parameter id="alpha" value="216.187851877211" constant="true"/>
+      <parameter id="beta" value="0.2" constant="false"/>
+      <parameter id="alpha0" value="0.216404256133345" constant="false"/>
+      <parameter id="alpha" value="216.187851877211" constant="false"/>
       <parameter id="eff" value="20" constant="true"/>
       <parameter id="n" value="2" constant="true"/>
       <parameter id="KM" value="40" constant="true"/>
       <parameter id="tau_mRNA" value="2" constant="true"/>
       <parameter id="tau_prot" value="10" constant="true"/>
-      <parameter id="t_ave" value="2.88539008177793" constant="true"/>
-      <parameter id="kd_mRNA" value="0.346573590279973" constant="true"/>
-      <parameter id="kd_prot" value="0.0693147180559945" constant="true"/>
-      <parameter id="k_tl" value="6.93147180559945" constant="true"/>
-      <parameter id="a_tr" value="29.97" constant="true"/>
+      <parameter id="t_ave" value="2.88539008177793" constant="false"/>
+      <parameter id="kd_mRNA" value="0.346573590279973" constant="false"/>
+      <parameter id="kd_prot" value="0.0693147180559945" constant="false"/>
+      <parameter id="k_tl" value="6.93147180559945" constant="false"/>
+      <parameter id="a_tr" value="29.97" constant="false"/>
       <parameter id="ps_a" value="0.5" constant="true"/>
       <parameter id="ps_0" value="0.0005" constant="true"/>
-      <parameter id="a0_tr" value="0.03" constant="true"/>
+      <parameter id="a0_tr" value="0.03" constant="false"/>
     </listOfParameters>
     <listOfRules>
       <assignmentRule variable="t_ave">
         <math xmlns="http://www.w3.org/1998/Math/MathML">
           <apply>
             <divide/>
-            <ci> _tau_mRNA </ci>
+            <ci> tau_mRNA </ci>
             <apply>
               <ln/>
               <cn type="integer"> 2 </cn>
@@ -54,17 +54,17 @@
         <math xmlns="http://www.w3.org/1998/Math/MathML">
           <apply>
             <divide/>
-            <ci> _tau_mRNA </ci>
-            <ci> _tau_prot </ci>
+            <ci> tau_mRNA </ci>
+            <ci> tau_prot </ci>
           </apply>
         </math>
       </assignmentRule>
       <assignmentRule variable="k_tl">
         <math xmlns="http://www.w3.org/1998/Math/MathML">
           <apply>
             <divide/>
-            <ci> _eff </ci>
-            <ci> _t_ave </ci>
+            <ci> eff </ci>
+            <ci> t_ave </ci>
           </apply>
         </math>
       </assignmentRule>
@@ -75,12 +75,12 @@
             <apply>
               <times/>
               <cn type="integer"> -60 </cn>
-              <ci> _ps_0 </ci>
+              <ci> ps_0 </ci>
             </apply>
             <apply>
               <times/>
               <cn type="integer"> 60 </cn>
-              <ci> _ps_a </ci>
+              <ci> ps_a </ci>
             </apply>
           </apply>
         </math>
@@ -90,7 +90,7 @@
           <apply>
             <times/>
             <cn type="integer"> 60 </cn>
-            <ci> _ps_0 </ci>
+            <ci> ps_0 </ci>
           </apply>
         </math>
       </assignmentRule>
@@ -102,7 +102,7 @@
               <ln/>
               <cn type="integer"> 2 </cn>
             </apply>
-            <ci> _tau_prot </ci>
+            <ci> tau_prot </ci>
           </apply>
         </math>
       </assignmentRule>
@@ -114,7 +114,7 @@
               <ln/>
               <cn type="integer"> 2 </cn>
             </apply>
-            <ci> _tau_mRNA </ci>
+            <ci> tau_mRNA </ci>
           </apply>
         </math>
       </assignmentRule>
@@ -124,13 +124,13 @@
             <divide/>
             <apply>
               <times/>
-              <ci> _a_tr </ci>
-              <ci> _eff </ci>
-              <ci> _tau_prot </ci>
+              <ci> a_tr </ci>
+              <ci> eff </ci>
+              <ci> tau_prot </ci>
             </apply>
             <apply>
               <times/>
-              <ci> _KM </ci>
+              <ci> KM </ci>
               <apply>
                 <ln/>
                 <cn type="integer"> 2 </cn>
@@ -145,13 +145,13 @@
             <divide/>
             <apply>
               <times/>
-              <ci> _a0_tr </ci>
-              <ci> _eff </ci>
-              <ci> _tau_prot </ci>
+              <ci> a0_tr </ci>
+              <ci> eff </ci>
+              <ci> tau_prot </ci>
             </apply>
             <apply>
               <times/>
-              <ci> _KM </ci>
+              <ci> KM </ci>
               <apply>
                 <ln/>
                 <cn type="integer"> 2 </cn>
@@ -220,6 +220,9 @@
         <listOfProducts>
           <speciesReference species="PX" stoichiometry="1" constant="false"/>
         </listOfProducts>
+        <listOfModifiers>
+          <modifierSpeciesReference species="X"/>
+        </listOfModifiers>
         <kineticLaw>
           <math xmlns="http://www.w3.org/1998/Math/MathML">
             <apply>
@@ -237,6 +240,9 @@
         <listOfProducts>
           <speciesReference species="PY" stoichiometry="1" constant="false"/>
         </listOfProducts>
+        <listOfModifiers>
+          <modifierSpeciesReference species="Y"/>
+        </listOfModifiers>
         <kineticLaw>
           <math xmlns="http://www.w3.org/1998/Math/MathML">
             <apply>
@@ -254,6 +260,9 @@
         <listOfProducts>
           <speciesReference species="PZ" stoichiometry="1" constant="false"/>
         </listOfProducts>
+        <listOfModifiers>
+          <modifierSpeciesReference species="Z"/>
+        </listOfModifiers>
         <kineticLaw>
           <math xmlns="http://www.w3.org/1998/Math/MathML">
             <apply>
@@ -322,6 +331,9 @@
         <listOfProducts>
           <speciesReference species="X" stoichiometry="1" constant="false"/>
         </listOfProducts>
+        <listOfModifiers>
+          <modifierSpeciesReference species="PZ"/>
+        </listOfModifiers>
         <kineticLaw>
           <math xmlns="http://www.w3.org/1998/Math/MathML">
             <apply>
@@ -363,6 +375,9 @@
         <listOfProducts>
           <speciesReference species="Y" stoichiometry="1" constant="false"/>
         </listOfProducts>
+        <listOfModifiers>
+          <modifierSpeciesReference species="PX"/>
+        </listOfModifiers>
         <kineticLaw>
           <math xmlns="http://www.w3.org/1998/Math/MathML">
             <apply>
@@ -404,6 +419,9 @@
         <listOfProducts>
           <speciesReference species="Z" stoichiometry="1" constant="false"/>
         </listOfProducts>
+        <listOfModifiers>
+          <modifierSpeciesReference species="PY"/>
+        </listOfModifiers>
         <kineticLaw>
           <math xmlns="http://www.w3.org/1998/Math/MathML">
             <apply>