Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fixed more SBML loading and writing bugs #25

Merged
merged 3 commits into from
Apr 23, 2020
Merged

Fixed more SBML loading and writing bugs #25

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
d08c64c
fixed multiple sbml bugs
ayush9pandey Apr 23, 2020
ee799b2
removed duplicate function getAllIds
ayush9pandey Apr 23, 2020
93200ad
minor bug fixed
ayush9pandey Apr 23, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
82 changes: 45 additions & 37 deletions bioscrape/sbmlutil.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand Down Expand Up @@ -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]
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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()
Expand All @@ -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")
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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
4 changes: 2 additions & 2 deletions bioscrape/types.pyx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand All @@ -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]
Expand Down
419 changes: 222 additions & 197 deletions examples/Basic Examples - START HERE.ipynb
View file
Open in desktop

Large diffs are not rendered by default.

72 changes: 45 additions & 27 deletions examples/models/represillator_sbml_back.xml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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>
Expand All @@ -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>
Expand All @@ -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>
Expand All @@ -90,7 +90,7 @@
<apply>
<times/>
<cn type="integer"> 60 </cn>
<ci> _ps_0 </ci>
<ci> ps_0 </ci>
</apply>
</math>
</assignmentRule>
Expand All @@ -102,7 +102,7 @@
<ln/>
<cn type="integer"> 2 </cn>
</apply>
<ci> _tau_prot </ci>
<ci> tau_prot </ci>
</apply>
</math>
</assignmentRule>
Expand All @@ -114,7 +114,7 @@
<ln/>
<cn type="integer"> 2 </cn>
</apply>
<ci> _tau_mRNA </ci>
<ci> tau_mRNA </ci>
</apply>
</math>
</assignmentRule>
Expand All @@ -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>
Expand All @@ -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>
Expand Down Expand Up @@ -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>
Expand All @@ -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>
Expand All @@ -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>
Expand Down Expand Up @@ -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>
Expand Down Expand Up @@ -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>
Expand Down Expand Up @@ -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>
Expand Down
Loading