# -*- coding: utf-8 -*-
"""
Created on Tue Dec 6 17:42:31 2011
@author: proto
"""
from copy import deepcopy, copy
from bionetgen.atomizer.writer import bnglWriter as writer
log = {"species": [], "reactions": []}
import re, sympy
from collections import Counter
from collections import defaultdict
import math as pymath
from bionetgen.atomizer.utils.util import logMess, TranslationException
import libsbml
from bionetgen.atomizer.bngModel import bngModel
import sympy
from sympy import Function
from sympy.abc import _clash
from sympy.printing.str import StrPrinter
from sympy.core.sympify import SympifyError
# Define 2 and 3 argument functions
# for sympy parsing
[docs]class sympyPiece(Function):
nargs = (3, 4, 5)
[docs]class sympyIF(Function):
nargs = 3
[docs]class sympyGT(Function):
nargs = 2
[docs]class sympyLT(Function):
nargs = 2
[docs]class sympyGEQ(Function):
nargs = 2
[docs]class sympyLEQ(Function):
nargs = 2
[docs]class sympyAnd(Function):
nargs = (2, 3, 4, 5)
[docs]class sympyOr(Function):
nargs = (2, 3, 4, 5)
[docs]class sympyNot(Function):
nargs = 1
[docs]def factorial(x):
temp = x
acc = 1
while temp > 0:
acc *= temp
temp -= 1
return acc
[docs]def comb(x, y, exact=True):
return factorial(x) / (factorial(y) * factorial(x - y))
bioqual = [
"BQB_IS",
"BQB_HAS_PART",
"BQB_IS_PART_OF",
"BQB_IS_VERSION_OF",
"BQB_HAS_VERSION",
"BQB_IS_HOMOLOG_TO",
"BQB_IS_DESCRIBED_BY",
"BQB_IS_ENCODED_BY",
"BQB_ENCODES",
"BQB_OCCURS_IN",
"BQB_HAS_PROPERTY",
"BQB_IS_PROPERTY_OF",
"BQB_HAS_TAXON",
"BQB_UNKNOWN",
]
modqual = [
"BQM_IS",
"BQM_IS_DESCRIBED_BY",
"BQM_IS_DERIVED_FROM",
"BQM_IS_INSTANCE_OF",
"BQM_HAS_INSTANCE",
"BQM_UNKNOWN",
]
annotationHeader = {"BQB": "bqbiol", "BQM": "bmbiol"}
[docs]def unrollSBMLFunction(function, sbmlFunctions):
"""
remove calls to functions inside functions
"""
modificationFlag = True
recursionIndex = 0
# remove calls to other sbml functions
while modificationFlag and recursionIndex < 20:
modificationFlag = False
for sbml in sbmlFunctions:
if sbml in function:
temp = writer.extendFunction(function, sbml, sbmlFunctions[sbml])
if temp != function:
function = temp
modificationFlag = True
recursionIndex += 1
break
return function
[docs]class SBML2BNGL:
"""
contains methods for extracting and formatting those sbml elements
that are translatable into bngl
"""
def __init__(self, model, useID=True, replaceLocParams=True, obs_map_file=None):
self.bngModel = bngModel()
self.bngModel.useID = useID
self.bngModel.replaceLocParams = replaceLocParams
self.bngModel.obs_map_file = obs_map_file
self.useID = useID
self.replaceLocParams = replaceLocParams
self.model = model
self.tags = {}
self.speciesUnits = {}
self.isConversion = True
self.boundaryConditionVariables = []
self.speciesDictionary = {}
self.speciesMemory = []
self.speciesAnnotationDict = None
self.reactionDictionary = {}
self.speciesAnnotation = None
self.speciesCompartments = None
self.unitDefinitions = self.getUnitDefinitions()
self.bngModel.unitDefinitions = self.unitDefinitions
self.convertSubstanceUnits = False
# keep track of obs names
self.obs_names = []
self.obs_map = {}
self.param_repl = {}
# ASS - I think there should be a check for compartments right here
# to determine if a) any compartment is actually used and
# b) remove unused ones and c) if only one is used with volume
# 1.0, remove that and write a normal model without compartments
# The following will give access to all methods to this information
self.noCompartment = None
# multi compartment warning flag
self.multi_comp_warn = False
self.check_noCompartment()
self.getSpecies()
# We want to re-use all the symbols
# generated by Sympy and then access it later on
self.all_syms = {}
self.all_syms.update(_clash)
self.all_syms.update({"pow": pow})
# self.all_syms["__epsilon__"] = sympy.symbols("__epsilon__")
# Do we want this? we'll have to handle it downstream
# handling a lot of functions via a dummy function
# that can take 1,2,3 variables
self.all_syms["Piecewise"] = sympyPiece
self.all_syms["piecewise"] = sympyPiece
self.all_syms["gt"] = sympyGT
self.all_syms["lt"] = sympyLT
self.all_syms["if"] = sympyIF
self.all_syms["sympyAnd"] = sympyAnd
self.all_syms["sympyOr"] = sympyOr
self.all_syms["sympyNot"] = sympyNot
# We are trying to replace things that we know
# are only in assignment rules in functions
self.only_assignment_dict = {}
if not hasattr(self, "used_molecules"):
self.used_molecules = []
if not hasattr(self, "used_symbols"):
self.used_symbols = []
if not hasattr(self, "arule_map"):
self.arule_map = {}
# Only write epsilon if we must
self.write_epsilon = False
[docs] def setConversion(self, conversion):
self.isConversion = conversion
[docs] def reset(self):
self.tags = {}
self.boundaryConditionVariables = []
self.speciesDictionary = {}
self.speciesMemory = []
self.getSpecies()
self.reactionDictionary = {}
[docs] def static_var(varname, value):
def decorate(func):
setattr(func, varname, value)
return func
return decorate
[docs] def isSameNameDifferentCompartment(self, name):
speciesList = []
for species in self.model.getListOfSpecies():
if species.getName() == name:
speciesList.append(species.getCompartment())
return len(speciesList) == len(set(speciesList))
[docs] def getRawSpecies(self, species, parameters=[], logEntries=True):
"""
*species* is the element whose SBML information we will extract
this method gets information directly
from an SBML related to a particular species.
It returns id,initialConcentration,(bool)isconstant and isboundary,
and the compartment
It also accounts for the fact that sometimes ppl use the same name for
molecules with different identifiers
"""
identifier = species.getId()
name = species.getName()
if name == "":
name = identifier
if species.isSetInitialConcentration():
initialValue = species.getInitialConcentration()
else:
initialValue = species.getInitialAmount()
isConstant = species.getConstant()
isBoundary = species.getBoundaryCondition()
# FIXME: this condition means that a variable/species can be changed
# by rules and/or events. this means that we effectively need a variable
# changed by a function that tracks this value, and all references
# to this observable have to be changed to the referrencing variable.
# http://sbml.org/Software/libSBML/docs/java-api/org/sbml/libsbml/Species.html
if isBoundary and not isConstant:
# isConstant = True
if (
not species.isSetInitialConcentration()
and not species.isSetInitialAmount()
):
initialValue = 1
compartment = species.getCompartment()
boundaryCondition = species.getBoundaryCondition()
standardizedName = standardizeName(name)
# if its a species that appends the compartment name remove it, it is not necessary in bionetgen
if standardizedName.endswith("_{0}".format(compartment)):
standardizedName = standardizedName.replace("_{0}".format(compartment), "")
# if standardizedName in ['Source','Trash','Sink']:
# standardizedName = '0'
if standardizedName in parameters:
standardizedName = "sp_{0}".format(standardizedName)
# it cannot start with a number
if standardizedName[:1].isdigit():
standardizedName = "s" + standardizedName
# two species cannot have the same name. Ids are unique but less informative, however typically species can be differentiated
# by compartment
if logEntries and standardizedName != "0":
if standardizedName in self.speciesMemory:
if len(list(self.model.getListOfCompartments())) == 1:
standardizedName += "_" + species.getId()
else:
# we can differentiate by compartment tag, no need to attach it to the name
# however an actual check needs to be made to make sure that the modeler is actually
# changing compartment information. If not, use id to differentiate.
if not self.isSameNameDifferentCompartment(species.getName()):
standardizedName += "_" + species.getId()
self.speciesMemory.append(standardizedName)
if boundaryCondition:
if standardizedName not in self.boundaryConditionVariables:
self.boundaryConditionVariables.append(standardizedName)
self.speciesDictionary[identifier] = standardizedName
returnID = identifier if self.useID else self.speciesDictionary[identifier]
if species.isSetSubstanceUnits():
self.speciesUnits[returnID] = species.getSubstanceUnits()
values = {}
values["returnID"] = returnID
if species.isSetInitialConcentration():
values["initialConcentration"] = initialValue
values["initialAmount"] = -1
elif species.isSetInitialAmount():
values["initialAmount"] = initialValue
values["initialConcentration"] = -1
else:
values["initialAmount"] = -1
values["initialConcentration"] = -1
values["isConstant"] = isConstant
values["isBoundary"] = isBoundary
values["compartment"] = compartment
values["name"] = name
values["identifier"] = identifier
if len(species.conversion_factor) > 0:
# we have a conversion factor to account for
values["conversionFactor"] = species.conversion_factor
else:
values["conversionFactor"] = None
return values
[docs] def getPrunnedTree(self, math, remainderPatterns, artificialObservables={}):
"""
walks through a series of * nodes and removes the remainder reactant factors
arg:remainderPatterns: argumetns to be removed from the tree
it also changes references to time variables to the keyword 'Time'
artificialObservables: species that are changed through an sbml assignment rule. their
usage in bng requires special handling.
"""
swapDict = {libsbml.AST_NAME_TIME: "Time"}
for node in [
x for x in (math.getLeftChild(), math.getRightChild()) if x != None
]:
if node.getType() in list(swapDict.keys()):
node.setName(swapDict[node.getType()])
if math.getCharacter() == "+" and math.getNumChildren() == 1:
math = math.getLeftChild()
while (math.getCharacter() == "*" or math.getCharacter() == "/") and len(
remainderPatterns
) > 0:
if libsbml.formulaToString(math.getLeftChild()) in remainderPatterns:
remainderPatterns.remove(libsbml.formulaToString(math.getLeftChild()))
if math.getCharacter() == "*":
math = math.getRightChild()
else:
math.getLeftChild().setValue(1)
elif libsbml.formulaToString(math.getRightChild()) in remainderPatterns:
remainderPatterns.remove(libsbml.formulaToString(math.getRightChild()))
math = math.getLeftChild()
else:
if (math.getLeftChild().getCharacter()) == "*":
math.replaceChild(
0, self.getPrunnedTree(math.getLeftChild(), remainderPatterns)
)
if (math.getRightChild().getCharacter()) == "*":
math.replaceChild(
math.getNumChildren() - 1,
self.getPrunnedTree(math.getRightChild(), remainderPatterns),
)
break
return math
[docs] def getIsTreeNegative(self, math):
"""
walks through a series of * nodes and detects whether there's a negative factor
fixme: we should actually test if the number of negative factors is odd
right now we are relying on the modelers not being malicious
when writing their rates laws.
"""
if math.getCharacter() == "*" or math.getCharacter() == "/":
if math.getLeftChild().isUMinus():
math.setCharacter("+")
# math.getLeftChild().setValue(long(libsbml.formulaToString(math.getLeftChild())[1:]))
return True
elif math.getLeftChild().getNumChildren() == 0 and libsbml.formulaToString(
math.getLeftChild()
).startswith("-"):
math.getLeftChild().setValue(
int(libsbml.formulaToString(math.getLeftChild())[1:])
)
return True
elif math.getRightChild().isUMinus():
math.setCharacter("+")
# math.getRightChild().setValue(long(libsbml.formulaToString(math.getRightChild())[1:]))
return True
elif math.getRightChild().getNumChildren() == 0 and libsbml.formulaToString(
math.getRightChild()
).startswith("-"):
math.getRightChild().setValue(
int(libsbml.formulaToString(math.getRightChild())[1:])
)
return True
else:
if (math.getLeftChild().getCharacter()) in ["*", "/", "-"]:
if self.getIsTreeNegative(math.getLeftChild()):
return True
if (math.getRightChild().getCharacter()) in ["*", "/", "-"]:
if self.getIsTreeNegative(math.getRightChild()):
return True
elif math.getCharacter() == "-" and math.getNumChildren() == 1:
math.setCharacter("+")
return True
return False
[docs] def getUnitDefinitions(self):
self.unitDictionary = {}
for unitDefinition in self.model.getListOfUnitDefinitions():
unitList = []
for unit in unitDefinition.getListOfUnits():
correctedUnit = libsbml.Unit.convertToSI(unit)
# unitList.append({'kind':unit.getKind(), 'scale':unit.getScale(),'multiplier':unit.getMultiplier(), 'exponent': unit.getExponent(), 'name':name})
for unit2 in correctedUnit.getListOfUnits():
name = (
unit.getName() if unit.getName() else unitDefinition.getName()
)
unitList.append(
{
"kind": unit2.getKind(),
"scale": unit2.getScale(),
"multiplier": unit2.getMultiplier(),
"exponent": unit2.getExponent(),
"name": name,
}
)
self.unitDictionary[unitDefinition.getId()] = unitList
return self.unitDictionary
[docs] def preProcessStoichiometry(self, reactants):
"""
checks for reactants with the same name in the reactant list. This
is mainly to account for reactants that dont have the stoichiometry
flag properly set and instead appear repeated
"""
uniqueReactantDict = defaultdict(int)
for reactant in reactants:
uniqueReactantDict[reactant[0]] += reactant[1]
return [(x, uniqueReactantDict[x]) for x in uniqueReactantDict]
[docs] def removeFactorFromMath(self, math, reactants, products, artificialObservables):
"""
it also adds symmetry factors. this checks for symmetry in the species names
artificialObservables: species names that are changed through assignment rules. their use requires special care when calculating a rate
"""
ifStack = Counter()
remainderPatterns = []
highStoichoiMetryFactor = 1
processedReactants = self.preProcessStoichiometry(reactants)
# ASS: I'm doing a hack, this is a flag to indicate
# that a species appears on both sides of a reaction
bothSides = False
for x in processedReactants:
# this is the symmtery factor for the rate constant
highStoichoiMetryFactor *= factorial(x[1])
y = [i[1] for i in products if i[0] == x[0]]
if len(y) > 0:
bothSides = True
y = y[0] if len(y) > 0 else 0
# TODO: check if this actually keeps the correct dynamics
# this is basically there to address the case where theres more products
# than reactants (synthesis)
if x[1] > y:
highStoichoiMetryFactor /= comb(int(x[1]), int(y), exact=True)
# print("HSMF comb: {}".format(highStoichoiMetryFactor))
for counter in range(0, int(x[1])):
remainderPatterns.append(x[0])
# for x in products:
# highStoichoiMetryFactor /= math.factorial(x[1])
# remainderPatterns = [x[0] for x in reactants]
# print remainderPatterns, artificialObservables.keys()
math = self.getPrunnedTree(math, remainderPatterns)
rateR = libsbml.formulaToString(math)
for element in remainderPatterns:
ifStack.update([element])
# This is here to denote if we are removing an extra factor
# of one of the elements (on top of the one that might or might
# not be there) so we have our element in the denominator
# hence the special treatment.
for element in ifStack:
if ifStack[element] > 1:
# ASS - removing if statement from functional rate definitions
# rateR = '{1}/({0}^{2} + __epsilon__)'.format(element, rateR, ifStack[element])
rateR = "{1}/({0}^{2})".format(element, rateR, ifStack[element])
self.write_epsilon = True
else:
# ASS - removing if statement from functional rate definitions
# rateR = '{1}/({0} + __epsilon__)'.format(element, rateR)
rateR = "{1}/({0})".format(element, rateR)
# self.write_epsilon = True
numFactors = max(math.getNumChildren(), len(ifStack))
if pymath.isinf(highStoichoiMetryFactor):
rateR = "{0} * 1e20".format(rateR)
logMess(
"ERROR:SIM204", 'Found usage of "inf" inside function {0}'.format(rateR)
)
elif highStoichoiMetryFactor != 1 and bothSides:
# ASS
# there is something wrong here, this multiplies regular
# rate constant by the highStoichiometry value and it's simply
# incorrect.
# Update: I think not multiplying is correct for most cases,
# I think this changes when we have a species on both sides
# of the reaction. Then, and only then, this parsing is relevant
# I believe.
# Update: hence the "bothSides" flag
rateR = "{0} * {1}".format(rateR, int(highStoichoiMetryFactor))
# we are adding a factor to the rate so we need to account for it when
# we are constructing the bngl equation (we dont want constrant expressions in there)
numFactors = max(1, numFactors)
# print("reactants and products: {} \n {} \n".format(reactants, products))
# print("HSMF: {}".format(highStoichoiMetryFactor))
# print("rateR: {}".format(rateR))
# print("numFactors: {}".format(numFactors))
return rateR, numFactors
[docs] def isAmount(self, reactantName):
for species in self.model.getListOfSpecies():
if species.getName() == reactantName:
if species.isSetInitialAmount():
return True
return False
[docs] def calculate_factor(self, react, prod, expr, removed):
"""
Calculates factors from reactants and products?
"""
ifStack = Counter()
remainderPatterns = []
highStoichoiMetryFactor = 1
processedReactants = self.preProcessStoichiometry(react)
# ASS: I'm doing a hack, this is a flag to indicate
# that a species appears on both sides of a reaction
bothSides = False
for x in processedReactants:
# this is the symmtery factor for the rate constant
highStoichoiMetryFactor *= factorial(x[1])
y = [i[1] for i in prod if i[0] == x[0]]
if len(y) > 0:
bothSides = True
y = y[0] if len(y) > 0 else 0
if x[1] > y:
highStoichoiMetryFactor /= comb(int(x[1]), int(y), exact=True)
for counter in range(0, int(x[1])):
remainderPatterns.append(x[0])
for rem in removed:
if rem in remainderPatterns:
remainderPatterns.remove(rem)
for element in remainderPatterns:
ifStack.update([element])
# Replace certain symbols s.t. libsbml.parseFormula gets it
expr = expr.replace("**", "^")
numFactors = max(libsbml.parseFormula(expr).getNumChildren(), len(ifStack))
if pymath.isinf(highStoichoiMetryFactor):
logMess(
"ERROR:SIM204",
'Found usage of "inf" inside function. This might be due to non-integer stoichiometry as well.',
)
elif highStoichoiMetryFactor != 1 and bothSides:
numFactors = max(1, numFactors)
return numFactors
[docs] def find_all_symbols(self, math, reactionID):
time_warn = False
curr_keys = list(self.all_syms.keys())
if math is None:
logMess(
"ERROR:SIM211",
f"Math for reaction with ID '{reactionID}' is not defined",
)
raise TranslationException(
f"ERROR:SIM211: Math for reaction with ID '{reactionID}' is not defined"
)
l = math.getListOfNodes()
replace_dict = {}
for inode in range(l.getSize()):
node = l.get(inode)
# Sympy doesn't like "def" in our string
name = node.getName()
if name == "def":
self.all_syms["def"] = sympy.symbols("__DEF__")
replace_dict["def"] = "__DEF__"
continue
if name == "lambda":
self.all_syms["lambda"] = sympy.symbols("__LAMBDA__")
replace_dict["lambda"] = "__LAMBDA__"
continue
if name == "as":
self.all_syms["as"] = sympy.symbols("__AS__")
replace_dict["as"] = "__AS__"
continue
if name == "del":
self.all_syms["del"] = sympy.symbols("__DEL__")
replace_dict["del"] = "__DEL__"
continue
if name == "time":
# Need to raise warning but do it
# once per reaction
if time_warn:
pass
else:
time_warn = True
if reactionID is not None:
logMess(
"WARNING:RATE004",
"Kinetic law of reaction {} is using the time function. The time functionality in BNG might not work as expected and throw off the translation.".format(
reactionID
),
)
else:
logMess(
"WARNING:RATE004",
"At least one reactions kinetic law is using the time function. The time functionality in BNG might not work as expected and throw off the translation.",
)
if name is not None:
if name not in curr_keys:
self.all_syms[name] = sympy.symbols(name)
# let's parse the formula and get non-numerical symbols
form = libsbml.formulaToString(math)
# If we need to replace anything
# TODO: Replace all of these with regexp
for it in replace_dict.items():
form = form.replace(it[0], it[1])
# Let's also pool this in used_symbols
for sym in self.all_syms.keys():
if sym not in self.used_symbols:
self.used_symbols.append(sym)
# Sympy doesn't allow and/not/or to be used
# outside what it deems to be acceptable
# TODO: Replace all of these with regexp
if "piecewise(" in form:
form = form.replace("piecewise(", "sympyPiece(")
replace_dict["piecewise"] = "sympyPiece"
if "gt(" in form:
form = form.replace("gt(", "sympyGT(")
replace_dict["gt"] = "sympyGT"
if "geq(" in form:
form = form.replace("geq(", "sympyGEQ(")
replace_dict["geq"] = "sympyGEQ"
if "lt(" in form:
form = form.replace("lt(", "sympyLT(")
replace_dict["lt"] = "sympyLT"
if "leq(" in form:
form = form.replace("leq(", "sympyLEQ(")
replace_dict["leq"] = "sympyLEQ"
if "if(" in form:
form = form.replace("if(", "sympyIF(")
replace_dict["if"] = "sympyIF"
if "and(" in form:
form = form.replace("and(", "sympyAnd(")
replace_dict["and"] = "sympyAnd"
# TODO: "or(" catches stuff like "floor(" and other
# potential functions. This needs to be extended
# to more potential or statements (e.g. *or(, +or( etc
# the same goes for other functions too but this is
# particularly a problem for this one
if " or(" in form:
form = form.replace("or(", "sympyOr(")
replace_dict["or"] = "sympyOr"
if "not(" in form:
form = form.replace("not(", "sympyNot(")
replace_dict["not"] = "sympyNot"
return form, replace_dict
[docs] def analyzeReactionRate(
self,
math,
compartmentList,
reversible,
rReactant,
rProduct,
reactionID,
parameterFunctions,
rModifier=[],
sbmlFunctions={},
split_rxn=False,
):
"""
This functions attempts to obtain the left and right hand sides of a rate reaction
function given a MathML tree. It also removes compartments and chemical factors from the function
Keyword arguments:
math -- the MathML math object
compartmentList -- a list of all the compartments in the system
reversible -- a boolean indicated whether there's should be a backward rate
rReactant -- a string list of the reactants.
rProduct -- a string list of the products.
sbmlFunctions -- a list of possible nested functiosn that we need to remove
"""
# We turn the AST object to string first
mathstring = libsbml.formulaToString(math)
# unroll sbml functions
mathstring = unrollSBMLFunction(mathstring, sbmlFunctions)
temp = libsbml.parseFormula(mathstring)
# if we could parse back after all those modifications...
if temp:
math = temp
# I'm removing the function defined and putting the code in here
# directly
react = rReactant
prod = rProduct
rateL = None
rateR = None
# Let's pull everything in the formula as symbols to use
# with sympify
form, replace_dict = self.find_all_symbols(math, reactionID)
# let's pull all names
all_names = [i[0] for i in react] + [i[0] for i in prod]
# SymPy is wonderful, _clash1 avoids built-ins like E, I etc
# FIXME:can we adjust the assignment rule stuff here?
try:
sym = sympy.sympify(form, locals=self.all_syms)
except SympifyError as e:
logMess(
"ERROR:SYMP001",
"Sympy couldn't parse a function, sorry but we can't handle this function definition.",
)
raise e # TranslationException(reactionID)
# Remove compartments if we use them.
# if not self.noCompartment:
compartments_to_remove = [sympy.symbols(comp) for comp in compartmentList]
# TODO: This is not fully correct, we need to know what
# compartment is on what side which is not currently
# being provided to this function
for comp in compartments_to_remove:
if comp in sym.atoms():
# Further issue, I know that this should be
# a multiplication but for BMD2 this is actually a
# problem? In fact, it looks like this is the case
# for regular mass action in SBML?
# This doesn't look right and it is a current
# hack?
n, d = sym.as_numer_denom()
if comp in n.atoms():
sym = sym / comp
elif comp in d.atoms():
sym = sym * comp
else:
pass
# If we are splitting, we don't need to do much
if split_rxn:
rate = str(sym).replace("**", "^")
for it in replace_dict.items():
rate = rate.replace(it[1], it[0])
return rate, "", 1, 1, False, split_rxn
# expand and take the terms out as left and right
exp = sympy.expand(sym)
# This shows if we can get X - Y
###### SPLIT RXN #######
# TODO: Figure out if something CAN be mass action
# and if not, just skip the rest and use split_rxn
###### SPLIT RXN #######
if exp.is_Add:
react_expr, prod_expr = self.gather_terms(exp)
# l,r = exp.as_two_terms()
# Let's also ensure that we have a + and - term
# if str(l).startswith("-") or str(r).startswith("-"):
# if str(l).startswith("-"):
# fwd_expr = r
# back_expr = l
# else:
# fwd_expr = l
# back_expr = r
if react_expr is None:
# this means we only have negative
# values in the kinetic laws. I guess
# this is possible if one of the parameters
# turns it into + overall rate. We'll turn this into
# unidirectional after multiplying by -, also will
# raise a warning
logMess(
"WARNING:RATE005",
"All terms in reaction {} are negative. It's possible that one of the parameters make the overall rate positive but ensure that this is what you expect".format(
reactionID
),
)
react_expr = -1 * prod_expr
prod_expr = None
if prod_expr is not None:
# Also get and parse the symbols
react_bols = [x[0] for x in react]
prod_bols = [x[0] for x in prod]
react_symbols = sympy.symbols(react_bols)
prod_symbols = sympy.symbols(prod_bols)
# Now we can manipulate it
# react_expr = fwd_expr
removedL = []
for ibol, bol in enumerate(react_symbols):
stoi = int(react[ibol][1])
# Now we can remove it
react_expr = react_expr / (bol**stoi)
removedL += [str(bol) for i in range(stoi)]
# Check if we can get 0 in the denominator
add_eps_react = False
n, d = react_expr.as_numer_denom()
for ibol, bol in enumerate(react_symbols):
if bol in d.atoms():
d = d.subs(bol, 0)
if d == 0:
# logMess('WARNING:RATE001', 'Denominator of rate constant in reaction {} can be 0. We are adding a small value epsilon to avoid discontinuities which can cause small errors in the model.'.format(reactionID))
add_eps_react = True
# let's instead split the rxn
split_rxn = True
rate = str(sym).replace("**", "^")
for it in replace_dict.items():
rate = rate.replace(it[1], it[0])
return rate, "", 1, 1, False, split_rxn
# prod_expr = back_expr
removedR = []
for ibol, bol in enumerate(prod_symbols):
stoi = int(prod[ibol][1])
# Now we can remove it
prod_expr = prod_expr / (bol**stoi)
removedR += [str(bol) for i in range(stoi)]
# Check if we can get 0 in the denominator
add_eps_prod = False
n, d = prod_expr.as_numer_denom()
for ibol, bol in enumerate(prod_symbols):
if bol in d.atoms():
d = d.subs(bol, 0)
if d == 0:
# logMess('WARNING:RATE001', 'Denominator of rate constant in reaction {} can be 0. We are adding a small value epsilon to avoid discontinuities which can cause small errors in the model.'.format(reactionID))
# let's instead split the rxn
split_rxn = True
rate = str(sym).replace("**", "^")
add_eps_prod = True
for it in replace_dict.items():
rate = rate.replace(it[1], it[0])
return rate, "", 1, 1, False, split_rxn
# prod_expr = prod_expr * -1
# TODO: We still need to figure out if we have
# our reactant/products in our expressions and
# if so set the nl/nr values accordingly
# Reproducing current behavior + expansion
re_proc = react_expr.nsimplify().evalf().simplify()
pe_proc = prod_expr.nsimplify().evalf().simplify()
# Adding epsilon if we have to
if add_eps_react:
# n,d = re_proc.as_numer_denom()
# rateL = "(" + str(n) + ")/(" + str(d) + "+__epsilon__)"
# self.write_epsilon = True
# add_eps_prod = True
# instead splitting the reaction
split_rxn = True
rate = str(sym).replace("**", "^")
for it in replace_dict.items():
rate = rate.replace(it[1], it[0])
return rate, "", 1, 1, False, split_rxn
else:
rateL = str(re_proc)
if add_eps_prod:
# n,d = pe_proc.as_numer_denom()
# rateR = "(" + str(n) + ")/(" + str(d) + "+__epsilon__)"
# self.write_epsilon = True
# add_eps_prod = True
# instead splitting the reaction
split_rxn = True
rate = str(sym).replace("**", "^")
for it in replace_dict.items():
rate = rate.replace(it[1], it[0])
return rate, "", 1, 1, False, split_rxn
else:
rateR = str(pe_proc)
nl = self.calculate_factor(react, prod, rateL, removedL)
nr = self.calculate_factor(prod, react, rateR, removedR)
# nl, nr = 2, 2
# BNG power function is ^ and not **
rateL = rateL.replace("**", "^")
rateR = rateR.replace("**", "^")
reversible = True
if rateL is None:
# if not simply reversible, rely on the SBML spec
if reversible:
logMess(
"WARNING:RATE002",
"SBML claims reversibility of reaction {} but the kinetic law is not easily separated to forward and backward rate constants. Assuming irreversible reaction instead.".format(
reactionID
),
)
# Also get and parse the symbols
react_bols = [x[0] for x in react]
react_symbols = sympy.symbols(react_bols)
# Now we can manipulate it
react_expr = exp
removedL = []
for ibol, bol in enumerate(react_symbols):
stoi = int(react[ibol][1])
# Now we can remove it
react_expr = react_expr / (bol**stoi)
removedL += [str(bol) for i in range(stoi)]
# Check if we can get the denominator to be 0
add_eps_react = False
n, d = react_expr.as_numer_denom()
for ibol, bol in enumerate(react_symbols):
if bol in d.atoms():
d = d.subs(bol, 0)
if d == 0:
# logMess('WARNING:RATE001', 'Denominator of rate constant in reaction {} can be 0. We are adding a small value epsilon to avoid discontinuities which can cause small errors in the model.'.format(reactionID))
# add_eps_react = True
# instead splitting the reaction
split_rxn = True
rate = str(sym).replace("**", "^")
for it in replace_dict.items():
rate = rate.replace(it[1], it[0])
return rate, "", 1, 1, False, split_rxn
re_proc = react_expr.nsimplify().evalf().simplify()
if add_eps_react:
# n,d = re_proc.as_numer_denom()
# rateL = "(" + str(n) + ")/(" + str(d) + "+__epsilon__)"
# self.write_epsilon = True
# instead splitting the reaction
split_rxn = True
rate = str(sym).replace("**", "^")
for it in replace_dict.items():
rate = rate.replace(it[1], it[0])
return rate, "", 1, 1, False, split_rxn
else:
rateL = str(re_proc)
nl = self.calculate_factor(react, prod, rateL, removedL)
rateL = rateL.replace("**", "^")
# Make unidirectional
rateR = "0"
nr = -1
reversible = False
# Check to ensure we don't have a negative
# rate constant
if rateL.startswith("-"):
logMess(
"WARNING:RATE003",
"After analysis, forward rate constant of reaction {} starts with - which is likely an error.".format(
reactionID
),
)
MrateL = rateL
MrateR = rateR
Mnl, Mnr = nl, nr
uRev = reversible
for it in replace_dict.items():
MrateL = MrateL.replace(it[1], it[0])
MrateR = MrateR.replace(it[1], it[0])
return MrateL, MrateR, Mnl, Mnr, uRev, split_rxn
def __getRawRules(
self, reaction, symmetryFactors, parameterFunctions, translator, sbmlfunctions
):
zerospecies = ["emptyset", "trash", "sink", "source"]
split_rxn = False
if self.useID:
reactant = [
(
reactant.getSpecies(),
reactant.getStoichiometry(),
reactant.getSpecies(),
)
for reactant in reaction.getListOfReactants()
if reactant.getSpecies().lower() not in zerospecies
and reactant.getStoichiometry() not in [0, "0"]
]
product = [
(product.getSpecies(), product.getStoichiometry(), product.getSpecies())
for product in reaction.getListOfProducts()
if product.getSpecies().lower() not in zerospecies
and product.getStoichiometry() not in [0, "0"]
]
else:
reactant = [
(
self.speciesDictionary[rElement.getSpecies()],
rElement.getStoichiometry(),
rElement.getSpecies(),
)
for rElement in reaction.getListOfReactants()
if self.speciesDictionary[rElement.getSpecies()].lower()
not in zerospecies
and rElement.getStoichiometry() not in [0, "0"]
]
product = [
(
self.speciesDictionary[rProduct.getSpecies()],
rProduct.getStoichiometry(),
rProduct.getSpecies(),
)
for rProduct in reaction.getListOfProducts()
if self.speciesDictionary[rProduct.getSpecies()].lower()
not in zerospecies
and rProduct.getStoichiometry() not in [0, "0"]
]
kineticLaw = reaction.getKineticLaw()
reversible = reaction.getReversible()
if kineticLaw is None:
return {
"reactants": reactant,
"products": product,
"parameters": [],
"rates": ["0", "0"],
"reversible": reversible,
"reactionID": reaction.getId(),
"numbers": [0, 0],
"modifiers": None,
"split_rxn": split_rxn,
}
rReactant = []
rProduct = []
# in case a given species was defined as the zero molecule don't include it in the rate correction method
for x in reaction.getListOfReactants():
# check if we are splitting the reaction due to conversion factors
if x.getSpecies() in self.bngModel.molecules:
if self.bngModel.molecules[x.getSpecies()].conversionFactor is not None:
split_rxn = True
if (
x.getSpecies().lower() not in zerospecies
and x.getStoichiometry()
not in [
0,
"0",
]
):
if not x.getConstant() and pymath.isnan(x.getStoichiometry()):
logMess(
"ERROR:SIM241",
"BioNetGen does not support non constant stoichiometries. Reaction {0} is not correctly translated".format(
reaction.getId()
),
)
raise TranslationException(
"ERROR:SIM241: BioNetGen does not support non constant stoichiometries. Reaction {0} is not correctly translated".format(
reaction.getId()
)
)
else:
speciesName = (
self.speciesDictionary[x.getSpecies()]
if x.getSpecies() in self.speciesDictionary
else ""
)
if (
speciesName in translator
and str(translator[speciesName]) == "0"
):
continue
rReactant.append((x.getSpecies(), x.getStoichiometry()))
for x in reaction.getListOfProducts():
# check if we are splitting the reaction due to conversion factors
if x.getSpecies() in self.bngModel.molecules:
if self.bngModel.molecules[x.getSpecies()].conversionFactor is not None:
split_rxn = True
if (
x.getSpecies().lower() not in zerospecies
and x.getStoichiometry()
not in [
0,
"0",
]
):
if not x.getConstant() and pymath.isnan(x.getStoichiometry()):
logMess(
"ERROR:SIM207",
"BioNetGen does not support non constant stoichiometries. Reaction {0} is not correctly translated".format(
reaction.getId()
),
)
raise TranslationException(
"ERROR:SIM207: BioNetGen does not support non constant stoichiometries. Reaction {0} is not correctly translated".format(
reaction.getId()
)
)
else:
speciesName = (
self.speciesDictionary[x.getSpecies()]
if x.getSpecies() in self.speciesDictionary
else ""
)
if (
speciesName in translator
and str(translator[speciesName]) == "0"
):
continue
rProduct.append((x.getSpecies(), x.getStoichiometry()))
# rReactant = [(x.getSpecies(), x.getStoichiometry()) for x in reaction.getListOfReactants() if x.getSpecies() not in ['EmptySet']]
# rProduct = [(x.getSpecies(), x.getStoichiometry()) for x in reaction.getListOfProducts() if x.getSpecies() not in ['EmptySet']]
rModifiers = [
x.getSpecies()
for x in reaction.getListOfModifiers()
if x.getSpecies() != "EmptySet"
]
parameters = [
(parameter.getId(), parameter.getValue(), parameter.getUnits())
for parameter in kineticLaw.getListOfParameters()
]
rateL = rateR = nl = nr = None
if True:
# TODO: For some reason creating a deepcopy of this screws everything up, even
# though its what we should be doing
# update: apparently the solution was to use copy instead of deepcopy. This is because
# the underlying swig code in c was causing conflicts when copied. make sure this actually works
math = copy(kineticLaw.getMath())
math = math.deepCopy()
# get a list of compartments so that we can remove them
compartmentList = []
for compartment in self.model.getListOfCompartments():
compartmentList.append(compartment.getId())
# check for non-integer stoichiometry
stoi = sum([r[1] for r in rReactant]) + sum([l[1] for l in rProduct])
if pymath.isnan(stoi):
logMess(
"ERROR:SIM208",
f"Stoichiometry for reaction {reaction.getId()} is NaN, can't translate.",
)
raise TranslationException(
f"ERROR:SIM241 Stoichiometry for reaction {reaction.getId()} is NaN, can't translate."
)
try:
if stoi != int(stoi):
non_integer_stoi = True
else:
non_integer_stoi = False
except:
logMess(
"ERROR:SIM209",
f"Stoichiometry for reaction {reaction.getId()} is {stoi}, can't translate.",
)
raise TranslationException(
f"ERROR:SIM208 Stoichiometry for reaction {reaction.getId()} is {stoi}, can't translate."
)
# remove compartments from expression. also separate left hand and right hand side
# check reaction molecularity and determine if it's
# possible to have a mass action to begin with
r_stoi = sum([r[1] for r in rReactant])
if r_stoi > 3 or non_integer_stoi:
# modified analyze reaction rate call because we'll be splitting
# the raection downstream
split_rxn = True
rateL, rateR, nl, nr, uReversible, split_rxn = self.analyzeReactionRate(
math,
compartmentList,
reversible,
rReactant,
rProduct,
reaction.getId(),
parameterFunctions,
rModifiers,
sbmlfunctions,
split_rxn,
)
if rateR == "0":
reversible = False
# FIXME: make sure this actually works
if symmetryFactors[0] > 1:
rateL = "({0})*({1})".format(rateL, symmetryFactors[0])
if symmetryFactors[1] > 1:
rateR = "({0})*({1})".format(rateR, symmetryFactors[1])
# we need to resolve observables BEFORE we do this
for obs_key in self.obs_map:
resL = re.search(r"(\W|^){0}(\W|$)".format(obs_key), rateL)
if resL is not None:
rateL = re.sub(
r"(\W|^)({0})(\W|$)".format(obs_key),
r"\g<1>{0}\g<3>".format(self.obs_map[obs_key]),
rateL,
)
resR = re.search(r"(\W|^){0}(\W|$)".format(obs_key), rateR)
if resR is not None:
rateR = re.sub(
r"(\W|^)({0})(\W|$)".format(obs_key),
r"\g<1>{0}\g<3>".format(self.obs_map[obs_key]),
rateR,
)
if not self.useID:
srateL = self.convertToName(rateL)
srateR = self.convertToName(rateR)
else:
srateL = rateL
srateR = rateR
if uReversible:
pass
# return compartments if the reaction is unimolecular
# they were removed in the first palce because its easier to handle
# around the equation in tree form when it has less terms
"""
if len(self.model.getListOfCompartments()) > 0:
for compartment in (self.model.getListOfCompartments()):
if compartment.getId() not in compartmentList:
if len(rReactant) != 2:
rateL = '{0} * {1}'.format(rateL,compartment.getSize())
if len(rProduct) != 2:
rateR = '{0} * {1}'.format(rateR,compartment.getSize())
"""
return {
"reactants": reactant,
"products": product,
"parameters": parameters,
"rates": [srateL, srateR],
"reversible": uReversible,
"reactionID": reaction.getId(),
"numbers": [nl, nr],
"modifiers": rModifiers,
"split_rxn": split_rxn,
}
[docs] def getReactionCenter(self, reactant, product, translator):
rcomponent = Counter()
pcomponent = Counter()
for element in reactant:
if element[0] in translator:
for molecule in translator[element[0]].molecules:
for component in molecule.components:
molecule.sort()
rcomponent.update(
Counter(
[
(
molecule.name,
component.name,
len(component.bonds) > 0,
component.activeState,
)
]
)
)
for element in product:
if element[0] in translator:
for molecule in translator[element[0]].molecules:
molecule.sort()
for component in molecule.components:
molecule.sort()
pcomponent.update(
Counter(
[
(
molecule.name,
component.name,
len(component.bonds) > 0,
component.activeState,
)
]
)
)
reactionCenter = [
(x[0], x[1])
for x in rcomponent
for y in pcomponent
if (x[0], x[1]) == (y[0], y[1])
and x != y
and rcomponent[x] != pcomponent[y]
]
rreactionCenter = [
(x[0], x[1])
for x in pcomponent
for y in rcomponent
if (x[0], x[1]) == (y[0], y[1])
and x != y
and pcomponent[x] != rcomponent[y]
]
return reactionCenter, rreactionCenter
[docs] def updateComponentCount(self, counterArray, reference, updateValue):
for element in counterArray:
if reference in counterArray[element]:
counterArray[element][reference] += updateValue
[docs] def reduceComponentSymmetryFactors(self, reaction, translator, functions):
"""
create symmetry factors for reactions with components and species with
identical names. This checks for symmetry in the components names then.
"""
# FIXME: This is entirely broken
zerospecies = ["emptyset", "trash", "sink", "source"]
if self.useID:
reactant = [
(rElement.getSpecies(), rElement.getStoichiometry())
for rElement in reaction.getListOfReactants()
if rElement.getSpecies() != "EmptySet"
]
product = [
(product.getSpecies(), product.getStoichiometry())
for product in reaction.getListOfProducts()
if product.getSpecies() != "EmptySet"
]
else:
reactant = [
(
self.speciesDictionary[rElement.getSpecies()],
rElement.getStoichiometry(),
)
for rElement in reaction.getListOfReactants()
]
product = [
(
self.speciesDictionary[rProduct.getSpecies()],
rProduct.getStoichiometry(),
)
for rProduct in reaction.getListOfProducts()
]
kineticLaw = reaction.getKineticLaw()
reversible = reaction.getReversible()
if kineticLaw is None:
return 1, 1
rReactant = rProduct = []
for x in reaction.getListOfReactants():
if (
x.getSpecies().lower() not in zerospecies
and x.getStoichiometry() not in [0, "0"]
and pymath.isnan(x.getStoichiometry())
):
if not x.getConstant():
logMess(
"ERROR:SIM241",
"BioNetGen does not support non constant stoichiometries. Reaction {0} is not correctly translated".format(
reaction.getId()
),
)
return 1, 1
else:
rReactant.append(x.getSpecies(), x.getStoichiometry())
for x in reaction.getListOfProducts():
if (
x.getSpecies().lower() not in zerospecies
and x.getStoichiometry() not in [0, "0"]
and pymath.isnan(x.getStoichiometry())
):
if not x.getConstant():
logMess(
"ERROR:SIM241",
"BioNetGen does not support non constant stoichiometries. Reaction {0} is not correctly translated".format(
reaction.getId()
),
)
return 1, 1
else:
rProduct.append(x.getSpecies(), x.getStoichiometry())
# TODO: For some reason creating a deepcopy of this screws everything up, even
# though its what we should be doing
rcomponent = defaultdict(Counter)
pcomponent = defaultdict(Counter)
# get the total count of components in the reactants and products
# e.g. components across diffent species
freactionCenter, breactionCenter = self.getReactionCenter(
reactant, product, translator
)
for element in reactant:
if element[0] in translator:
for molecule in translator[element[0]].molecules:
for component in molecule.components:
molecule.sort()
componentList = Counter([(molecule.signature(freactionCenter))])
for _ in range(0, int(element[1])):
rcomponent[
(
molecule.name,
component.name,
len(component.bonds) > 0,
component.activeState,
)
].update(componentList)
for element in product:
if element[0] in translator:
for molecule in translator[element[0]].molecules:
molecule.sort()
for component in molecule.components:
componentList = Counter([(molecule.signature(breactionCenter))])
for _ in range(0, int(element[1])):
pcomponent[
(
molecule.name,
component.name,
len(component.bonds) > 0,
component.activeState,
)
].update(componentList)
"""
only keep information for reaction centers
"""
reactionCenters = [
(x[0], x[1])
for x in rcomponent
for y in pcomponent
if (x[0], x[1]) == (y[0], y[1]) and x != y
]
rcomponent = {
x: rcomponent[x] for x in rcomponent if (x[0], x[1]) in reactionCenters
}
pcomponent = {
x: pcomponent[x] for x in pcomponent if (x[0], x[1]) in reactionCenters
}
# is the number of components across products and reactants the same?
# eg is there any DeleteMolecules action
pcorrectionFactor = 1
rcorrectionFactor = 1
rStack = []
pStack = []
"""
if a reaction can take place in several ways account for it in the reaction
rate (this is specially important in dimer and trimer binding)
pcomponent[element] < rcomponent[element] asks if an specific instance
of a component decreases in number from a reactant to a product
for example if there are 3 A(b)'s and one binds, we will have 2 A(b)'s
in the product
"""
rcomponentTemp = deepcopy(rcomponent)
pcomponentTemp = deepcopy(pcomponent)
# calculate actual symmetry factors
for key in rcomponent:
if key in pcomponent:
for element in rcomponent[key]:
if rcomponent[key] == 1:
continue
# if theres a component on one side of the equation that
# appears a different number of times on the other side of the equation
if element in pcomponent[key]:
if (
pcomponent[key][element] < rcomponent[key][element]
and set([key[0].lower(), key[1].lower()]) not in rStack
):
rcorrectionFactor *= comb(
rcomponent[key][element],
pcomponent[key][element],
exact=1,
)
rStack.append(set([key[0].lower(), key[1].lower()]))
# once we choose a component for a previous action
# this limits the options for subsequent actions
# although substracting one from the target sites
# may not be the right option. double check.
self.updateComponentCount(pcomponent, element, -1)
else:
for element2 in pcomponent[key]:
if (
pcomponent[key][element2] < rcomponent[key][element]
and set([key[0].lower(), key[1].lower()]) not in rStack
):
rcorrectionFactor *= comb(
rcomponent[key][element],
pcomponent[key][element2],
exact=1,
)
rStack.append(set([key[0].lower(), key[1].lower()]))
self.updateComponentCount(pcomponent, element2, -1)
rcomponent = rcomponentTemp
pcomponent = pcomponentTemp
if reversible:
for key in pcomponent:
if key in rcomponent:
for element in pcomponent[key]:
if pcomponent[key] == 1:
continue
if element in rcomponent[key]:
if (
rcomponent[key][element] < pcomponent[key][element]
and set([key[0].lower(), key[1].lower()]) not in pStack
):
pcorrectionFactor *= comb(
pcomponent[key][element],
rcomponent[key][element],
exact=1,
)
pStack.append(set([key[0].lower(), key[1].lower()]))
self.updateComponentCount(rcomponent, element, -1)
else:
for element2 in rcomponent[key]:
if (
rcomponent[key][element2] < pcomponent[key][element]
and set([key[0].lower(), key[1].lower()])
not in pStack
):
pcorrectionFactor *= comb(
pcomponent[key][element],
rcomponent[key][element2],
exact=1,
)
pStack.append(set([key[0].lower(), key[1].lower()]))
self.updateComponentCount(rcomponent, element2, -1)
return rcorrectionFactor, pcorrectionFactor
[docs] def convertToName(self, rate):
for element in sorted(self.speciesDictionary, key=len, reverse=True):
while True:
oldRate = rate
if element in rate:
rate = re.sub(
r"(\W|^)({0})(\W|$)".format(element),
r"\1{0}\3".format(self.speciesDictionary[element]),
rate,
)
if rate != oldRate:
continue
break
# rate = rate.replace(element,self.speciesDictionary[element])
return rate
def __getRawCompartments(self, compartment):
"""
Private method used by the getCompartments method
"""
idid = compartment.getId()
name = compartment.getName()
size = compartment.getSize()
# volume messes up the reactions
# size = 1.0
dimensions = compartment.getSpatialDimensions()
if dimensions in [0, 1]:
logMess(
"WARNING:SIM103",
"{1}-D compartments are not supported. Changing for 2-D compartments for {0}. Please verify this does not affect simulation".format(
name, dimensions
),
)
dimensions = 2
if not compartment.getConstant():
# this compartment is not constant in size, not supported right now
logMess(
"ERROR:SIM203",
f"Compartments with variable size is not supported, compartment: {idid}",
)
return idid, dimensions, size, name
def __getRawFunctions(self, function):
math = function[1].getMath()
name = function[1].getId()
return name, libsbml.formulaToString(math)
[docs] def getSBMLFunctions(self):
functions = {}
for function in enumerate(self.model.getListOfFunctionDefinitions()):
functionInfo = self.__getRawFunctions(function)
if functionInfo[1] is None:
logMess(
"ERROR:SIM212",
f"Math for function with ID '{functionInfo[0]}' is not defined",
)
raise TranslationException(
f"ERROR:SIM212: Math for function with ID '{functionInfo[0]}' is not defined"
)
functions[functionInfo[0]] = writer.bnglFunction(
functionInfo[1],
functionInfo[0],
[],
reactionDict=self.reactionDictionary,
)
return functions
[docs] def getCompartments(self):
"""
Returns an array of triples, where each triple is defined as
(compartmentName,dimensions,size)
"""
compartments = []
# ASS - return empty if we are not using compartments
units = None
if self.noCompartment:
return compartments
unitDefinitions = self.getUnitDefinitions()
if "volume" in unitDefinitions:
units = "#volume units: {0}".format(
"*".join([x["name"] for x in unitDefinitions["volume"]])
)
compartments.append(units)
else:
units = "#volume units: L"
compartments.append(units)
for _, compartment in enumerate(self.model.getListOfCompartments()):
compartmentInfo = self.__getRawCompartments(compartment)
comp_obj = self.bngModel.make_compartment()
# ASS - removing "cell" as default compartment
name = compartmentInfo[0]
# making the obj before adding
comp_obj.Id = name
comp_obj.dim = compartmentInfo[1]
comp_obj.size = compartmentInfo[2]
comp_obj.unit = units
if name != compartmentInfo[3] and compartmentInfo[3] != "":
compartments.append(
"%s %d %s #%s"
% (name, compartmentInfo[1], compartmentInfo[2], compartmentInfo[3])
)
comp_obj.cmt = compartmentInfo[3]
else:
compartments.append(
"%s %d %s" % (name, compartmentInfo[1], compartmentInfo[2])
)
self.bngModel.add_compartment(comp_obj)
return compartments
[docs] def updateFunctionReference(self, reaction, updatedReferences):
newRate = reaction[3]
for reference in updatedReferences:
newRate = re.sub(
r"(\W|^)({0})(\W|$)".format(reference),
r"\1{0}\3".format(updatedReferences[reference]),
newRate,
)
return newRate
[docs] def getSymmetryFactors(self, reaction):
zerospecies = ["emptyset", "trash", "sink", "source"]
if self.useID:
reactant = [
(rElement.getSpecies(), rElement.getStoichiometry())
for rElement in reaction.getListOfReactants()
if rElement.getSpecies() != "EmptySet"
]
product = [
(product.getSpecies(), product.getStoichiometry())
for product in reaction.getListOfProducts()
if product.getSpecies() != "EmptySet"
]
else:
reactant = [
(
self.speciesDictionary[rElement.getSpecies()],
rElement.getStoichiometry(),
)
for rElement in reaction.getListOfReactants()
]
product = [
(
self.speciesDictionary[rProduct.getSpecies()],
rProduct.getStoichiometry(),
)
for rProduct in reaction.getListOfProducts()
]
react_counts = {}
for react in reactant:
if react[0] in react_counts:
react_counts[react[0]] += react[1]
else:
react_counts[react[0]] = react[1]
if len(react_counts) == 0:
lfact = 1
else:
lfact = max(react_counts.values())
prod_counts = {}
for prod in product:
if prod[0] in prod_counts:
prod_counts[prod[0]] += prod[1]
else:
prod_counts[prod[0]] = prod[1]
if len(prod_counts) == 0:
rfact = 1
else:
rfact = max(prod_counts.values())
return lfact, rfact
[docs] def getReactions(
self,
translator={},
isCompartments=False,
extraParameters={},
atomize=False,
parameterFunctions={},
database=None,
):
"""
@returns: a triple containing the parameters,reactions,functions
"""
# @FIXME:this part of the code is there so that we only generate the functions list once through different
# iterations of this call. This is because we cannot create a clone of the 'math' object for this
# reaction and it is being permanently changed every call. It's ugly but it works. Change for something
# better when we figure out how to clone the math object
if not hasattr(self.getReactions, "functionFlag"):
self.getReactions.__func__.functionFlag = False or (not atomize)
reactions = []
reactionStructure = []
parameters = []
functions = []
# We want to keep track of the molecules/species we
# actually used in the reactions
functionTitle = "fRate"
# self.unitDefinitions = self.getUnitDefinitions()
database.rawreactions = []
if len(self.model.getListOfReactions()) == 0:
logMess(
"WARNING:SIM104",
"Model contains no natural reactions, all reactions are produced by SBML rules",
)
# At the end of this loop, we want to go in and adjust any cases
# where a value that can go to 0 is in the denominator
for index, reaction in enumerate(self.model.getListOfReactions()):
parameterDict = {}
currParamConv = {}
# symmetry factors for components with the same name
# FIXME: This reduceComponentSymmetryFactors is completely broken
# and will only give 1,1 right now
# sl, sr = self.reduceComponentSymmetryFactors(
# reaction, translator, functions
# )
sl, sr = self.getSymmetryFactors(reaction)
sbmlfunctions = self.getSBMLFunctions()
try:
rawRules = self.__getRawRules(
reaction, [sl, sr], parameterFunctions, translator, sbmlfunctions
)
database.rawreactions.append(rawRules)
except TranslationException as e:
if database and database.ignore:
reactions.append(
"#Reaction {0} is not correctly translated. check log for details".format(
e.value
)
)
continue
else:
raise TranslationException(e.value + " during reaction processing")
# making the rule object
rule_obj = self.bngModel.make_rule()
rule_obj.rule_ind = index + 1
rule_obj.parse_raw(rawRules)
rule_obj.symm_factors = [sl, sr]
# Let's add our molecules
for r in rawRules["reactants"]:
if r[0] not in self.used_molecules:
self.used_molecules.append(r[0])
for p in rawRules["products"]:
if p[0] not in self.used_molecules:
self.used_molecules.append(p[0])
if len(rawRules["parameters"]) > 0:
for parameter in rawRules["parameters"]:
"""
if self.convertSubstanceUnits:
#newParameter = self.convertToStandardUnits(parameter[1],self.unitDictionary)
parameters.append('r%d_%s %f' % (index+1, parameter[0], parameter[1]))
parameterDict[parameter[0]] = parameter[1]
else:
"""
curr_param_name = "r%d_%s" % (index + 1, parameter[0])
curr_param = "r%d_%s %f" % (index + 1, parameter[0], parameter[1])
parameters.append(curr_param)
parameterDict[parameter[0]] = parameter[1]
currParamConv[parameter[0]] = curr_param_name
# add these to parameters in obj
param_obj = self.bngModel.make_parameter()
param_obj.Id = curr_param_name
param_obj.val = parameter[1]
param_obj.rxn_ind = index + 1
self.bngModel.add_parameter(param_obj)
compartmentList = []
compartmentList.extend(
[
[self.__getRawCompartments(x)[0], self.__getRawCompartments(x)[2]]
for x in self.model.getListOfCompartments()
]
)
threshold = 0
if rule_obj.raw_num[0] > threshold or rule_obj.raw_rates[0] in translator:
functionName = "%s%d()" % (functionTitle, index)
else:
# append reactionNumbers to parameterNames
finalString = str(rule_obj.raw_rates[0])
for parameter in parameterDict:
finalString = re.sub(
r"(\W|^)({0})(\W|$)".format(parameter),
r"\1{0}\3".format("r{0}_{1}".format(index + 1, parameter)),
finalString,
)
functionName = finalString
if self.getReactions.functionFlag and "delay" in rule_obj.raw_rates[0]:
logMess(
"ERROR:SIM202",
"BNG cannot handle delay functions in function %s" % functionName,
)
fobj = self.bngModel.make_function()
fobj.Id = functionName
fobj.rule_ptr = rule_obj
fobj.compartmentList = compartmentList
if rule_obj.reversible:
if (
rule_obj.raw_num[0] > threshold
or rule_obj.raw_rates[0] in translator
):
fobj.definition = rule_obj.raw_rates[0]
if self.getReactions.functionFlag:
# local parameter replacement flag
if self.replaceLocParams:
fstr = writer.bnglFunction(
rule_obj.raw_rates[0],
functionName,
rawRules["reactants"],
compartmentList,
parameterDict,
self.reactionDictionary,
)
functions.append(fstr)
fobj.local_dict = parameterDict
self.bngModel.add_function(fobj)
else:
fstr = writer.bnglFunction(
rule_obj.raw_rates[0],
functionName,
rawRules["reactants"],
compartmentList,
currParamConv,
self.reactionDictionary,
)
functions.append(fstr)
fobj.local_dict = currParamConv
self.bngModel.add_function(fobj)
if (
rawRules["numbers"][1] > threshold
or rule_obj.raw_rates[1] in translator
):
functionName2 = "%s%dm()" % (functionTitle, index)
fobj_2 = self.bngModel.make_function()
fobj_2.Id = functionName2
fobj_2.rule_ptr = rule_obj
fobj_2.definition = rule_obj.raw_rates[1]
fobj_2.compartmentList = compartmentList
if self.getReactions.functionFlag:
# local parameter replacement flag
if self.replaceLocParams:
functions.append(
writer.bnglFunction(
rule_obj.raw_rates[1],
functionName2,
rule_obj.raw_prod,
compartmentList,
parameterDict,
self.reactionDictionary,
)
)
fobj_2.local_dict = parameterDict
self.bngModel.add_function(fobj_2)
else:
functions.append(
writer.bnglFunction(
rule_obj.raw_rates[1],
functionName2,
rule_obj.raw_prod,
compartmentList,
currParamConv,
self.reactionDictionary,
)
)
fobj_2.local_dict = currParamConv
self.bngModel.add_function(fobj_2)
self.reactionDictionary[
rawRules["reactionID"]
] = "({0} - {1})".format(functionName, functionName2)
finalRateStr = "{0},{1}".format(functionName, functionName2)
rule_obj.rate_cts = (functionName, functionName2)
else:
finalString = str(rawRules["rates"][1])
for parameter in parameterDict:
finalString = re.sub(
r"(\W|^)({0})(\W|$)".format(parameter),
r"\1{0}\3".format("r{0}_{1}".format(index + 1, parameter)),
finalString,
)
finalRateStr = "{0},{1}".format(functionName, finalString)
rule_obj.rate_cts = (functionName, finalString)
else:
if (
rawRules["numbers"][0] > threshold
or rawRules["rates"][0] in translator
):
fobj.definition = rule_obj.raw_rates[0]
if self.getReactions.functionFlag:
# local parameter replacement flag
if self.replaceLocParams:
functions.append(
writer.bnglFunction(
rawRules["rates"][0],
functionName,
rawRules["reactants"],
compartmentList,
parameterDict,
self.reactionDictionary,
)
)
fobj.local_dict = parameterDict
self.bngModel.add_function(fobj)
else:
functions.append(
writer.bnglFunction(
rawRules["rates"][0],
functionName,
rawRules["reactants"],
compartmentList,
currParamConv,
self.reactionDictionary,
)
)
fobj.local_dict = currParamConv
self.bngModel.add_function(fobj)
self.reactionDictionary[rawRules["reactionID"]] = "{0}".format(
functionName
)
finalRateStr = functionName
rule_obj.rate_cts = (functionName,)
reactants = [x for x in rawRules["reactants"]]
products = [x for x in rawRules["products"]]
modifierComment = (
"#Modifiers({0})".format(", ".join(rawRules["modifiers"]))
if rawRules["modifiers"]
else ""
)
#### ADD RXN SEP HERE ####
if rule_obj.raw_splt:
# if we are splitting, we want to remove symmetry factor
# from a fRate because it will spread over all reactants
# and products
if "fRate" in rule_obj.rate_cts[0]:
# we have functional rate constant
if int(rule_obj.symm_factors[0]) != 1:
# we have a non-zero symmetry factor
if rule_obj.rate_cts[0] not in self.bngModel.functions:
logMess(
"ERROR:SIM206",
"Rate constant function needs adjusting but can't find function: {}"
% functionName,
)
defn = self.bngModel.functions[rule_obj.rate_cts[0]].definition
self.bngModel.functions[
rule_obj.rate_cts[0]
].definition = f"({defn})/({rule_obj.symm_factors[0]})"
if rule_obj.reversible:
logMess(
"ERROR:SIM205",
"Splitting a reversible reaction, please check if correct, function: {}"
% functionName,
)
ctr = 0
# Now we write a single reaction for each
# member with modified reaction rate constants
# first RHS
for reactant in reactants:
r = reactant
stoi = r[1]
if int(stoi) != 1.0:
nRateStr = "-1*{}*({})".format(stoi, finalRateStr)
else:
nRateStr = "-1*({})".format(finalRateStr)
nr = (r[0], 1.0, r[2])
# adjust reaction name
rxn_name = rawRules["reactionID"] + "_reactants_" + str(ctr)
rxn_str = writer.bnglReaction(
[],
[nr],
nRateStr,
self.tags,
translator,
(
isCompartments
or (
(len(reactants) == 0 or len(products) == 0)
and self.getReactions.__func__.functionFlag
)
),
rawRules["reversible"],
reactionName=rxn_name,
comment=modifierComment,
)
reactions.append(rxn_str)
# same thing for the model
nrule_obj = self.bngModel.make_rule()
nrule_obj.parse_raw(rawRules)
nrule_obj.reversible = False
nrule_obj.Id = rxn_name
nrule_obj.rate_cts = (nRateStr,)
nrule_obj.reactants = []
nrule_obj.products = [nr]
self.bngModel.add_rule(nrule_obj)
# tick the ctr
ctr += 1
# then LHS
ctr = 0
for product in products:
p = product
stoi = p[1]
if int(stoi) != 1.0:
nRateStr = "{}*{}".format(stoi, finalRateStr)
else:
nRateStr = "{}".format(finalRateStr)
np = (p[0], 1.0, p[2])
# adjust reaction name
rxn_name = rawRules["reactionID"] + "_products_" + str(ctr)
rxn_str = writer.bnglReaction(
[],
[np],
nRateStr,
self.tags,
translator,
(
isCompartments
or (
(len(reactants) == 0 or len(products) == 0)
and self.getReactions.__func__.functionFlag
)
),
rawRules["reversible"],
reactionName=rxn_name,
comment=modifierComment,
)
reactions.append(rxn_str)
# same thing for the model
nrule_obj = self.bngModel.make_rule()
nrule_obj.parse_raw(rawRules)
nrule_obj.reversible = False
nrule_obj.Id = rxn_name
nrule_obj.rate_cts = (nRateStr,)
nrule_obj.reactants = []
nrule_obj.products = [np]
self.bngModel.add_rule(nrule_obj)
# tick the ctr
ctr += 1
#### END RXN SEP ####
else:
# add the rule
rule_obj.reactants = reactants
rule_obj.products = products
self.bngModel.add_rule(rule_obj)
rxn_str = writer.bnglReaction(
reactants,
products,
finalRateStr,
self.tags,
translator,
(
isCompartments
or (
(len(reactants) == 0 or len(products) == 0)
and self.getReactions.__func__.functionFlag
)
),
rawRules["reversible"],
reactionName=rawRules["reactionID"],
comment=modifierComment,
)
reactions.append(rxn_str)
if atomize:
self.getReactions.__func__.functionFlag = True
self.bngModel.tags = self.tags
return parameters, reactions, functions
[docs] def gather_terms(self, exp):
pos, neg = [], []
l, r = exp.as_two_terms()
resolve = [l, r]
while len(resolve) > 0:
elem = resolve.pop()
if elem.is_Add:
l, r = elem.as_two_terms()
resolve += [l, r]
else:
# TODO: Do we have a better check?
if str(elem).startswith("-"):
neg.append(elem)
else:
pos.append(elem)
# FIXME: Return None correctly
l, r = None, None
if len(pos) > 0:
l = pos.pop(0)
if len(pos) > 0:
for e in pos:
l += e
if len(neg) > 0:
r = -1 * neg.pop(0)
if len(neg) > 0:
for e in neg:
r += -1 * e
return l, r
def __getRawAssignmentRules(self, arule):
"""
Function to pull out rate constants (if going to be encoded
as a reaction) and properties of an assignment rule to be
parsed downstream.
"""
variable = arule.getVariable()
var = sympy.symbols(variable)
if variable not in self.all_syms.keys():
self.all_syms[variable] = var
# Sympy stuff
form, replace_dict = self.find_all_symbols(arule.getMath(), None)
# We might need this for debugging complicated
# piecewise function forms.
try:
sym = sympy.sympify(form, locals=self.all_syms)
except SympifyError as e:
logMess(
"ERROR:SYMP001",
"Sympy couldn't parse a function, sorry but we can't handle this function definition.",
)
raise TranslationException(arule.getId())
# if it's an assignment, it's going to be
# encoded as a unidirectional rxn
if not arule.isAssignment():
# expand and take the terms out as left and right
exp = sympy.expand(sym)
rateL = None
rateR = None
# This shows if we can get X - Y
if exp.is_Add:
react_expr, prod_expr = self.gather_terms(exp)
if react_expr is None:
# TODO: LogMess this
print("no forward reaction rate?")
# Let's also ensure that we have a + and - term
elif prod_expr is not None:
# Remove mass action
prod_expr = prod_expr / var
# Check for epsilon
add_eps_prod = False
n, d = prod_expr.as_numer_denom()
for bol in d.atoms():
d = d.subs(bol, 0)
if d == 0:
add_eps_prod = True
# Reproducing current behavior + expansion
re_proc = react_expr.nsimplify().evalf().simplify()
pe_proc = prod_expr.nsimplify().evalf().simplify()
# Adding epsilon if we have to
rateL = str(re_proc)
if add_eps_prod:
rateL = str(sym)
rateR = "0"
else:
rateR = str(pe_proc)
rateL = rateL.replace("**", "^")
rateR = rateR.replace("**", "^")
# slap back in the replacement stuff
for it in replace_dict.items():
rateL = rateL.replace(it[1], it[0])
rateR = rateR.replace(it[1], it[0])
if rateL is None:
# if not simply reversible, rely on the SBML spec
react_expr = exp
re_proc = react_expr.nsimplify().evalf().simplify()
rateL = str(re_proc)
rateL = rateL.replace("**", "^")
# Make unidirectional
rateR = "0"
# slap back in the replacement stuff
for it in replace_dict.items():
rateL = rateL.replace(it[1], it[0])
else:
rateL = str(sym)
rateR = "0"
# slap back in the replacement stuff
for it in replace_dict.items():
rateL = rateL.replace(it[1], it[0])
if not self.useID:
rateL = self.convertToName(rateL)
rateR = self.convertToName(rateR)
return variable, [rateL, rateR], arule.isAssignment(), arule.isRate()
[docs] def adjustInitialConditions(
self,
parameters,
initialConditions,
artificialObservables,
observables,
functions,
):
"""
assignment rules require further adjustment after parsed
to their initial values.
"""
# We'll need this to write nicely
prnter = StrPrinter({"full_prec": False})
# This gets the initial conditions we have to adjust
initCondsToAdjust = [
x for x in artificialObservables.keys() if x.endswith("_ar")
]
# return original values if nothing needs adjusted
if len(initCondsToAdjust) == 0:
return initialConditions
logMess("DEBUG:ARUL001", "adjusting assignment rule initial conditions")
# These are the formulas we need to calculate dynamically
formulaToAdjustWith = {}
for cond in initCondsToAdjust:
form = artificialObservables[cond]
name = form.split("=")[0].split()[0]
for func in functions:
splt = func.split("=")
if len(splt) > 2:
logMess(
"WARNING:ARUL002",
"There are multiple = signs in function {}".format(splt[0]),
)
if name == splt[0].split()[0]:
# This avoids issues with multiple = signs
formulaToAdjustWith[name] = splt[-1]
# Now we have the formulas to dynamically calculate, let's get a parameter dictionary to replace parameters down the line
param_dict = dict([(x.split()[0], x.split()[1]) for x in parameters])
# we need to parse observables to map
# names to species descriptions
obs_map = {}
for iobs, obs in enumerate(observables):
splt = obs.split()
obs_map[splt[2]] = splt[1]
# Now that we have observables, we can map the names
# to their initial values to be replaced later
initValMap = {}
initCondSplit = []
for initCond in initialConditions:
splt = initCond.split()
initCondSplit.append(splt)
# I'm a bit vary of this, not sure if this is
# the only way the $ might appear honestly
# keep an eye out for bugs here
if splt[0].startswith("$"):
check_name = splt[0][1:]
else:
check_name = splt[0]
# if the name is in the observable species defs
if check_name in obs_map.keys():
# we slap that into our initial value map
initValMap[obs_map[check_name]] = splt[1]
# This is to allow us to just use .split to be able to
# separate species/parameter stuff from math stuff
sympy_forms = {}
for item in formulaToAdjustWith.items():
name, form = item
try:
sympy_forms[name] = sympy.sympify(form, locals=self.all_syms)
except:
sympy_forms[name] = None
sympy_subs = {}
# We need to replace the values in sympy formulas
for item in sympy_forms.items():
name, form = item
if form is not None:
# Looping over parameters and replacing
for par_item in param_dict.items():
pn, pv = par_item
pv = pv.replace("(", "").replace(")", "")
try:
ps, pvs = sympy.symbols(pn), sympy.Number(pv)
except:
logMess(
"ERROR:SYMP003",
"Sympy couldn't parse an initial condition: {} and {}, this might cause issues for the initial state of your model".format(
pn, pv
),
)
continue
form = form.subs(ps, pvs)
# Replacing species from initial conditions
for spec_item in initValMap.items():
sn, sv = spec_item
sv = sv.replace("(", "").replace(")", "")
try:
ss, svs = sympy.symbols(sn), sympy.Number(sv)
except:
logMess(
"ERROR:SYMP003",
"Sympy couldn't parse an initial condition: {} and {}, this might cause issues for the initial state of your model".format(
sn, sv
),
)
continue
form = form.subs(ss, svs)
# And now converting the rest to zeros so we have a value
# that we can fully evaluate
for atm in form.atoms(sympy.Symbol):
form = form.subs(atm, 0)
sympy_subs[name] = form.evalf()
else:
sympy_subs[name] = 0
# Evaluate them using Sympy
adjustedInitialValues = []
for item in sympy_subs.items():
name, val = item
try:
val_str = prnter.doprint(val)
if val_str == "nan":
val = 0
adjustedInitialValues.append(val)
except:
adjustedInitialValues.append(0)
# Re-write initial conditions
adjustedInitialConditions = []
for iic, initCond in enumerate(initCondSplit):
toAdjustCheck = initCond[-1].replace("#", "") + "_ar"
if toAdjustCheck in initCondsToAdjust:
val = adjustedInitialValues[initCondsToAdjust.index(toAdjustCheck)]
initCondSplit[iic][1] = "{}".format(prnter.doprint(val))
adjustedInitialConditions.append(" ".join(initCond))
else:
adjustedInitialConditions.append(" ".join(initCond))
return adjustedInitialConditions
[docs] def getAssignmentRules(
self, zparams, parameters, molecules, observablesDict, translator
):
"""
this method obtains an SBML rate rules and assignment rules. They
require special handling since rules are often both defined as rules
and parameters initialized as 0, so they need to be removed from the parameters list
"""
# FIXME: This function removes compartment info and this leads to mis-replacement of variables downstream. e.g. Calc@ER and Calc@MIT both gets written as Calc and downstream the replacement is wrong.
# FIXME: This function gets a list of observables which sometimes are turned into assignment rules but then are not updated in the observablesDict. E.g. X_comp1 gets in, X_ar is created and you can't have BOTH X_comp1 in a reaction AND X_ar adjusting X itself. You MUST pick one, if both are happening raise and error and exit out. For now I'll say if we have _ar then we replace the X_comp1 with X_ar and test.
# Going to use this to match names and remove params
# if need be
param_map = dict([(x.split()[0], x) for x in parameters])
compartmentList = []
compartmentList.extend(
[
[self.__getRawCompartments(x)[0], self.__getRawCompartments(x)[2]]
for x in self.model.getListOfCompartments()
]
)
arules = []
aParameters = {}
zRules = zparams
removeParameters = []
artificialReactions = []
artificialObservables = {}
nonamecounter = 0
for arule in self.model.getListOfRules():
rawArule = self.__getRawAssignmentRules(arule)
# print("looping over assignment rules")
arule_obj = self.bngModel.make_arule()
arule_obj.parse_raw(rawArule)
arule_obj.compartmentList = compartmentList
# rule has no name
if arule_obj.Id == "":
logMess(
"ERROR:SIM215",
"atomizer has found an sbml rule without a name. {0}".format(
rawArule[1:]
),
)
rawArule = list(rawArule)
new_arule_name = "noname{0}".format(nonamecounter)
rawArule[0] = new_arule_name
nonamecounter += 1
arule_obj.Id = new_arule_name
if arule_obj.isRate == True:
# it is a rate rule
if rawArule[0] in self.boundaryConditionVariables:
logMess(
"WARNING:SIM105",
"rate rules ({0}) are not properly supported in BioNetGen simulator".format(
rawArule[0]
),
)
rateLaw1 = arule_obj.rates[0]
rateLaw2 = arule_obj.rates[1]
# TODO: Add to bngModel functions
arate_name = "arRate{0}".format(rawArule[0])
func_str = writer.bnglFunction(
rateLaw1,
arate_name,
[],
compartments=compartmentList,
reactionDict=self.reactionDictionary,
)
arules.append(func_str)
if rateLaw2 != "0":
# TODO: Add to bngModel functions
armrate_name = "armRate{0}".format(rawArule[0])
func2_str = writer.bnglFunction(
rateLaw2,
armrate_name,
[],
compartments=compartmentList,
reactionDict=self.reactionDictionary,
)
arules.append(func2_str)
# ASS2019 - I'm not sure if this is the right place to fix the tags. Basically, up until this point, the artificial reactions don't have tags. This results in the 0 <-> A type reactions to lack a compartment, leading to a non-functional BNGL file. I think the better solution might be during rule (SBML rule, not BNGL rule) parsing and update the parser/SBML2BNGL tags instead.
try:
comp = self.tags[rawArule[0]]
except KeyError:
# ASS - We need to default to an existing compartment if we are going
# to remove @cell as a default compartment
if not self.noCompartment:
self.tags[rawArule[0]] = "@" + compartmentList[0][0]
# ASS - If self.useID is set, use the ID value, not the name
if self.useID:
self.used_molecules.append(rawArule[0])
if rateLaw2 == "0":
rxn_str = writer.bnglReaction(
[],
[[rawArule[0], 1, rawArule[0]]],
"{0}".format("arRate{0}".format(rawArule[0])),
self.tags,
translator,
isCompartments=True,
comment="#rateLaw",
reversible=False,
)
else:
rxn_str = writer.bnglReaction(
[],
[[rawArule[0], 1, rawArule[0]]],
"{0},{1}".format(
"arRate{0}".format(rawArule[0]),
"armRate{0}".format(rawArule[0]),
),
self.tags,
translator,
isCompartments=True,
comment="#rateLaw",
)
artificialReactions.append(rxn_str)
else:
self.used_molecules.append(self.convertToName(rawArule[0]).strip())
if rateLaw2 == "0":
rxn_str = writer.bnglReaction(
[],
[[self.convertToName(rawArule[0]).strip(), 1, rawArule[0]]],
"{0}".format("arRate{0}".format(rawArule[0])),
self.tags,
translator,
isCompartments=True,
comment="#rateLaw",
reversible=False,
)
else:
rxn_str = writer.bnglReaction(
[],
[[self.convertToName(rawArule[0]).strip(), 1, rawArule[0]]],
"{0},{1}".format(
"arRate{0}".format(rawArule[0]),
"armRate{0}".format(rawArule[0]),
),
self.tags,
translator,
isCompartments=True,
comment="#rateLaw",
)
artificialReactions.append(rxn_str)
if rawArule[0] in zparams:
removeParameters.append("{0} 0".format(rawArule[0]))
zRules.remove(rawArule[0])
else:
for element in parameters:
# TODO: if for whatever reason a rate rule
# was defined as a parameter that is not 0
# remove it. This might not be exact behavior
if re.search("^{0}\s".format(rawArule[0]), element):
logMess(
"WARNING:SIM106",
"Parameter {0} corresponds both as a non zero parameter \
and a rate rule, verify behavior".format(
element
),
)
removeParameters.append(element)
# it is an assigment rule
elif arule_obj.isAssignment:
"""
normal species observables references in functions keep the format <speciesName>_<compartment> in function references,
and observables dict keeps track of that. however when a species is defined by an assignment function we wish to
keep track of reference <speciesName> that points to a standard BNGL function
"""
# it was originially defined as a zero parameter, so delete it from the parameter list definition
if rawArule[0] in zRules:
# dont show assignment rules as parameters
zRules.remove(rawArule[0])
matches = [
molecules[x]
for x in molecules
if molecules[x]["name"] == rawArule[0]
]
if matches:
if matches[0]["isBoundary"]:
artificialObservables[
rawArule[0] + "_ar"
] = writer.bnglFunction(
rawArule[1][0],
rawArule[0] + "_ar()",
[],
compartments=compartmentList,
reactionDict=self.reactionDictionary,
)
self.arule_map[rawArule[0]] = rawArule[0] + "_ar"
if rawArule[0] in observablesDict:
observablesDict[rawArule[0]] = rawArule[0] + "_ar"
continue
else:
logMess(
"ERROR:SIM201",
"Variables that are both changed by an assignment rule and reactions are not \
supported in BioNetGen simulator. The variable will be split into two".format(
rawArule[0]
),
)
artificialObservables[
rawArule[0] + "_ar"
] = writer.bnglFunction(
rawArule[1][0],
rawArule[0] + "_ar()",
[],
compartments=compartmentList,
reactionDict=self.reactionDictionary,
)
self.arule_map[rawArule[0]] = rawArule[0] + "_ar"
if rawArule[0] in observablesDict:
observablesDict[rawArule[0]] = rawArule[0] + "_ar"
continue
elif rawArule[0] in [observablesDict[x] for x in observablesDict]:
artificialObservables[
rawArule[0] + "_ar"
] = writer.bnglFunction(
rawArule[1][0],
rawArule[0] + "_ar()",
[],
compartments=compartmentList,
reactionDict=self.reactionDictionary,
)
self.arule_map[rawArule[0]] = rawArule[0] + "_ar"
if rawArule[0] in observablesDict:
observablesDict[rawArule[0]] = rawArule[0] + "_ar"
continue
elif rawArule[0] in molecules:
if molecules[rawArule[0]]["isBoundary"]:
# We should probably re-write this with the name since that's what's used other places
name = molecules[rawArule[0]]["returnID"]
artificialObservables[name + "_ar"] = writer.bnglFunction(
rawArule[1][0],
name + "_ar()",
[],
compartments=compartmentList,
reactionDict=self.reactionDictionary,
)
self.arule_map[rawArule[0]] = name + "_ar"
# TODO: Let's store what we know are assignment rules. We can maybe assume that, if something has an assignment rule, it can't in turn be in a reaction? If this is wrong, we can't model this anyway, so we should probably just make an assumption and let people know.
self.only_assignment_dict[name] = name + "_ar"
self.bngModel.add_arule(arule_obj)
continue
else:
# if not boundary but is a species, Jose
# is turning this into an assignment rule
# with a different name (uses ID).
# It looks as if the goal was to handle
# both situations via renaming.
# FIXME: This is very likely broken but
# I'm not 100% sure how it breaks things.
# TODO: Check, if we have this in observables we need to adjust the observablesDict because we are writing an assignment rule for this instead
name = molecules[rawArule[0]]["returnID"]
artificialObservables[name + "_ar"] = writer.bnglFunction(
rawArule[1][0],
name + "_ar()",
[],
compartments=compartmentList,
reactionDict=self.reactionDictionary,
)
self.arule_map[rawArule[0]] = name + "_ar"
self.only_assignment_dict[name] = name + "_ar"
if name in observablesDict:
observablesDict[name] = name + "_ar"
self.bngModel.add_arule(arule_obj)
continue
else:
# check if it is defined as an observable
# FIXME: This doesn't check for parameter namespace
# TODO: What is going on here?
candidates = [
idx for idx, x in enumerate(observablesDict) if rawArule[0] == x
]
assigObsFlag = False
for idx in candidates:
# if re.search('\s{0}\s'.format(rawArule[0]),observables[idx]):
artificialObservables[
rawArule[0] + "_ar"
] = writer.bnglFunction(
rawArule[1][0],
rawArule[0] + "_ar()",
[],
compartments=compartmentList,
reactionDict=self.reactionDictionary,
)
self.arule_map[rawArule[0]] = rawArule[0] + "_ar"
assigObsFlag = True
break
if assigObsFlag:
continue
# if its not a param/species/observable
# TODO: now, if we replace this with the returnID do we
# overlap with it's other possible uses?
# name = molecules[rawArule[0]]['returnID']
# self.only_assignment_dict[name] = name+"_ar"
# artificialObservables[name+'_ar'] = writer.bnglFunction(rawArule[1][0],name+'()',[],compartments=compartmentList,reactionDict=self.reactionDictionary)
# This doesn't actually check for clashes with
# parameter namespace
if rawArule[0] in param_map.keys():
removeParameters.append(param_map[rawArule[0]])
artificialObservables[rawArule[0] + "_ar"] = writer.bnglFunction(
rawArule[1][0],
rawArule[0] + "()",
[],
compartments=compartmentList,
reactionDict=self.reactionDictionary,
)
self.arule_map[rawArule[0]] = rawArule[0] + "_ar"
else:
"""
if for whatever reason you have a rule that is not assigment
or rate and it is initialized as a non zero parameter, give it
a new name
"""
if rawArule[0] not in zparams:
ruleName = "ar" + rawArule[0]
else:
ruleName = rawArule[0]
zRules.remove(rawArule[0])
arules.append(
writer.bnglFunction(
rawArule[1][0],
ruleName,
[],
compartments=compartmentList,
reactionDict=self.reactionDictionary,
)
)
aParameters[rawArule[0]] = "ar" + rawArule[0]
"""
elif rawArule[2] == True:
for parameter in parameters:
if re.search('^{0}\s'.format(rawArule[0]),parameter):
print '////',rawArule[0]
"""
# we can't decide any of this here, we need the
# full model setup before deciding
self.bngModel.add_arule(arule_obj)
return (
aParameters,
arules,
zRules,
artificialReactions,
removeParameters,
artificialObservables,
)
[docs] def convertToStandardUnits(self, parameterValue, unitDefinition):
for factor in unitDefinition:
if factor["scale"] != 0:
parameterValue *= 10 ** factor["scale"]
if factor["multiplier"] != 1:
parameterValue *= factor["multiplier"]
if factor["exponent"] != 1:
parameterValue **= factor["exponent"]
return parameterValue
[docs] def convertToStandardUnitString(self, parameterValue, unitDefinition):
for factor in unitDefinition:
if factor["multiplier"] != 1:
parameterValue = "({0} * {1})".format(
parameterValue, factor["multiplier"]
)
if factor["exponent"] != 1:
parameterValue = "({0} ^ {1})".format(
parameterValue, factor["exponent"]
)
if factor["scale"] != 0:
parameterValue = "({0} * 1e{1})".format(parameterValue, factor["scale"])
return parameterValue
def __getRawParameters(self, parameter):
parameterSpecs = {}
parameterSpecs["id"] = parameter.getId()
parameterSpecs["value"] = parameter.getValue()
parameterSpecs["name"] = parameter.getName()
parameterSpecs["units"] = parameter.getUnits()
return parameterSpecs
[docs] def getParameters(self):
parameters = []
zparam = []
for parameter in self.model.getListOfParameters():
parameterSpecs = (
parameter.getId(),
parameter.getValue(),
parameter.getConstant(),
parameter.getUnits(),
)
# reserved keywords
param_obj = self.bngModel.make_parameter()
if parameterSpecs[0] == "e":
# TODO: raise a warning
parameterSpecs = ("__e__", parameterSpecs[1])
self.param_repl["e"] = "__e__"
if parameterSpecs[1] == 0:
zparam.append(parameterSpecs[0])
else:
if parameter.getUnits() != "":
parameters.append(
"{0} {1} #units:{2}".format(
parameterSpecs[0], parameterSpecs[1], parameter.getUnits()
)
)
param_obj.units = parameter.getUnits()
else:
parameters.append(
"{0} {1}".format(parameterSpecs[0], parameterSpecs[1])
)
if len(parameterSpecs) == 3:
param_obj.Id = parameterSpecs[0]
param_obj.val = parameterSpecs[1]
param_obj.cts = parameterSpecs[2]
else:
param_obj.Id = parameterSpecs[0]
param_obj.val = parameterSpecs[1]
param_obj.cts = True
self.bngModel.add_parameter(param_obj)
self.bngModel.param_repl = self.param_repl
return parameters, zparam
# ASS - trying to implement a more complete solution for
# non compartmental models
[docs] def check_noCompartment(self, parameters=[]):
self.compartmentDict = {}
# self.compartmentDict[""] = 1
comps = [i for i in self.model.getListOfCompartments()]
if len(comps) == 0:
self.noCompartment = True
self.bngModel.noCompartment = True
return
for compartment in self.model.getListOfCompartments():
self.compartmentDict[compartment.getId()] = compartment.getSize()
self.noCompartment = False
self.bngModel.noCompartment = False
# Get all rawSpecies
allRawSpecies = map(
lambda x: self.getRawSpecies(x, parameters), self.model.getListOfSpecies()
)
# Pull all used compartments in those speices and make a set
allUsedCompartments = set(map(lambda x: x["compartment"], allRawSpecies))
# Basically, if we are using a single compartment of volume 1,
# we are not actually using a compartment, not from the POV of
# BNGL. To improve readability, we want to just write a regular
# BNGL model instead of a cBNGL model. Especially true since
# this is the case for most SBML models.
if len(allUsedCompartments) == 1:
# We are using only 1 compartment, check volume
# FIXME: We will try removing the compartment
# if only one is used
# self.noCompartment = True
# self.bngModel.noCompartment = True
if self.compartmentDict[allUsedCompartments.pop()] == 1:
# we have 1 compartment and it's volume is 1
# just don't use compartments.
self.noCompartment = True
self.bngModel.noCompartment = True
elif len(allUsedCompartments) > 1:
if not self.multi_comp_warn:
logMess(
"WARNING:COMP001",
"Multiple compartments are used, please note that Atomizer does not automatically try to infer your compartment topology which is important for how rules fire in cBNGL. Make sure your comparment topology is set correctly after translation",
)
self.multi_comp_warn = True
self.speciesMemory = []
[docs] def getSpecies(self, translator={}, parameters=[]):
"""
in sbml parameters and species have their own namespace. not so in
bionetgen, so we need to rename things if they share the same name
"""
def default_to_regular(d):
if isinstance(d, defaultdict):
d = {k: default_to_regular(v) for k, v in d.items()}
return d
# gotta reset the bngModel everytime
# this function is called
self.bngModel._reset()
# find concentration units
unitDefinitions = self.unitDefinitions
if "substance" in unitDefinitions:
substance = "*".join([x["name"] for x in unitDefinitions["substance"]])
else:
substance = "mol"
if "volume" in unitDefinitions:
volume = "/".join([x["name"] for x in unitDefinitions["volume"]])
else:
volume = "L"
concentrationUnits = "{0}/{1}".format(substance, volume)
annotationMoleculesText = {}
moleculesText = []
speciesText = []
observablesText = []
observablesDict = {}
names = []
rawSpeciesName = list(translator.keys())
speciesTranslationDict = {}
compartmentDict = {}
compartmentDict[""] = 1
speciesAnnotationInfo = default_to_regular(self.getFullAnnotation())
annotationInfo = {"moleculeTypes": {}, "species": {}}
for compartment in self.model.getListOfCompartments():
compartmentDict[compartment.getId()] = compartment.getSize()
unitFlag = True
for species in self.model.getListOfSpecies():
# making molecule and seed species objs for
# the obj based model
molec_obj = self.bngModel.make_molecule()
spec_obj = self.bngModel.make_species()
#
rawSpecies = self.getRawSpecies(species, parameters)
# letting the objs parse the rawSpecies
molec_obj.parse_raw(rawSpecies)
spec_obj.parse_raw(rawSpecies)
if rawSpecies["compartment"] != "":
# ASS - First change for "noCompartments"
if self.noCompartment:
rawSpecies["compartment"] = ""
self.tags[rawSpecies["identifier"]] = ""
else:
self.tags[rawSpecies["identifier"]] = "@%s" % (
rawSpecies["compartment"]
)
if rawSpecies["returnID"] in translator:
if rawSpecies["returnID"] in rawSpeciesName:
rawSpeciesName.remove(rawSpecies["returnID"])
if (
translator[rawSpecies["returnID"]].getSize() == 1
and translator[rawSpecies["returnID"]].molecules[0].name
not in names
and translator[rawSpecies["returnID"]].molecules[0].name
not in rawSpeciesName
):
names.append(translator[rawSpecies["returnID"]].molecules[0].name)
annotationTemp = []
if rawSpecies["returnID"] in speciesAnnotationInfo:
for annotation in speciesAnnotationInfo[rawSpecies["returnID"]]:
parts = annotation.split("_")
header = annotationHeader[parts[0]]
qual = parts[1].lower() + "".join(
[x.capitalize() for x in parts[2:]]
)
entry = ", ".join(
[
":".join(x.split("/")[-2:])
for x in speciesAnnotationInfo[
rawSpecies["returnID"]
][annotation]
]
)
annotationTemp.append(
"#^ {0}:{1} {2}".format(header, qual, entry)
)
# we'll add this to our model
self.bngModel.add_molecule(molec_obj)
mtext = translator[rawSpecies["returnID"]].str2()
moleculesText.append(mtext)
if rawSpecies["returnID"] in speciesAnnotationInfo:
annotationInfo["moleculeTypes"][
translator[rawSpecies["returnID"]].str2()
] = annotationTemp
del speciesAnnotationInfo[rawSpecies["returnID"]]
# TODO: Not sure if there are more examples of this
# but glucose in 380 has both a normal species AND
# a boundary species separately
# elif rawSpecies['isBoundary']:
# self.bngModel.add_molecule(molec_obj)
else:
# we'll add this to our model
self.bngModel.add_molecule(molec_obj)
mtext = rawSpecies["returnID"] + "()"
moleculesText.append(mtext)
if rawSpecies["returnID"] in speciesAnnotationInfo:
annotationInfo["moleculeTypes"][
rawSpecies["returnID"]
] = speciesAnnotationInfo[rawSpecies["returnID"]]
del speciesAnnotationInfo[rawSpecies["returnID"]]
# if rawSpecies['identifier'] == 'glx' and len(translator) > 0:
temp = "$" if rawSpecies["isConstant"] != 0 else ""
tmp = (
translator[str(rawSpecies["returnID"])]
if rawSpecies["returnID"] in translator
else rawSpecies["returnID"] + "()"
)
# this determines the name to be written
if (
rawSpecies["initialConcentration"] > 0
or rawSpecies["initialAmount"] > 0
):
tmp2 = temp
if rawSpecies["identifier"] in self.tags:
tmp2 = self.tags[rawSpecies["identifier"]]
if rawSpecies["initialAmount"] > 0.0:
# Removing the compartment section if we are not using it
if self.noCompartment:
speciesText.append(
"{1}{2} {3} #{4} #{5}".format(
tmp2,
temp,
str(tmp),
rawSpecies["initialAmount"],
rawSpecies["returnID"],
rawSpecies["identifier"],
)
)
else:
speciesText.append(
"{0}:{1}{2} {3} #{4} #{5}".format(
tmp2,
temp,
str(tmp),
rawSpecies["initialAmount"],
rawSpecies["returnID"],
rawSpecies["identifier"],
)
)
elif rawSpecies["initialConcentration"] > 0.0:
if self.isConversion:
# convert to molecule counts
if "substance" in unitDefinitions:
newParameterStr = self.convertToStandardUnitString(
rawSpecies["initialConcentration"],
unitDefinitions["substance"],
)
newParameter = self.convertToStandardUnits(
rawSpecies["initialConcentration"],
unitDefinitions["substance"],
) # conversion to moles
else:
newParameter = rawSpecies["initialConcentration"]
newParameterStr = str(rawSpecies["initialConcentration"])
newParameter = (
newParameter * 6.022e23
) # convertion to molecule counts
# get compartment size
if self.noCompartment:
compartmentSize = 1.0
else:
compartmentSize = self.model.getCompartment(
rawSpecies["compartment"]
).getSize()
newParameter = compartmentSize * newParameter
# temp testing AS
spec_obj.val = newParameter
spec_obj.isConc = False
# temp testing AS
if unitFlag:
if self.noCompartment:
speciesText.append(
"{1}{2} {3} # {4}mol/L * 6.022e23/mol *{7}L #{5} #{6}".format(
tmp2,
temp,
str(tmp),
newParameter,
newParameterStr,
rawSpecies["returnID"],
rawSpecies["identifier"],
compartmentSize,
concentrationUnits,
)
)
else:
speciesText.append(
"{0}:{1}{2} {3} # {4}mol/L * 6.022e23/mol *{7}L #{5} #{6}".format(
tmp2,
temp,
str(tmp),
newParameter,
newParameterStr,
rawSpecies["returnID"],
rawSpecies["identifier"],
compartmentSize,
concentrationUnits,
)
)
unitFlag = False
else:
if self.noCompartment:
speciesText.append(
"{1}{2} {3} #original {4}{8} #{5} #{6}".format(
tmp2,
temp,
str(tmp),
newParameter,
rawSpecies["initialConcentration"],
rawSpecies["returnID"],
rawSpecies["identifier"],
compartmentSize,
concentrationUnits,
)
)
else:
speciesText.append(
"{0}:{1}{2} {3} #original {4}{8} #{5} #{6}".format(
tmp2,
temp,
str(tmp),
newParameter,
rawSpecies["initialConcentration"],
rawSpecies["returnID"],
rawSpecies["identifier"],
compartmentSize,
concentrationUnits,
)
)
else:
if self.noCompartment:
speciesText.append(
"{1}{2} {3} #{4} #{5}".format(
tmp2,
temp,
str(tmp),
rawSpecies["initialConcentration"],
rawSpecies["returnID"],
rawSpecies["identifier"],
)
)
else:
speciesText.append(
"{0}:{1}{2} {3} #{4} #{5}".format(
tmp2,
temp,
str(tmp),
rawSpecies["initialConcentration"],
rawSpecies["returnID"],
rawSpecies["identifier"],
)
)
elif rawSpecies["isConstant"]:
if self.noCompartment:
speciesText.append(
"{1}{2} {3} #{4} #{5}".format(
tmp2,
temp,
str(tmp),
0,
rawSpecies["returnID"],
rawSpecies["identifier"],
)
)
else:
speciesText.append(
"{0}:{1}{2} {3} #{4} #{5}".format(
tmp2,
temp,
str(tmp),
0,
rawSpecies["returnID"],
rawSpecies["identifier"],
)
)
self.bngModel.add_species(spec_obj)
if rawSpecies["returnID"] == "e":
modifiedName = "__e__"
else:
modifiedName = rawSpecies["returnID"]
# user defined zero molecules are not included in the observable list
if str(tmp) != "0":
if (
rawSpecies["compartment"] != ""
and len(list(self.model.getListOfCompartments())) > 1
):
self.obs_names.append(modifiedName)
# self.obs_map[rawSpecies["identifier"]] = "{0}_{1}".format(
# modifiedName, rawSpecies["compartment"]
# )
# observablesText.append(
# "Species {0}_{3} @{3}:{1} #{2}".format(
# modifiedName,
# tmp,
# rawSpecies["name"],
# rawSpecies["compartment"],
# )
# )
# observablesDict[modifiedName] = "{0}_{1}".format(
# modifiedName, rawSpecies["compartment"]
# )
self.obs_map[rawSpecies["identifier"]] = "{0}".format(
modifiedName, rawSpecies["compartment"]
)
observablesText.append(
"Species {0} @{3}:{1} #{2}".format(
modifiedName,
tmp,
rawSpecies["name"],
rawSpecies["compartment"],
)
)
observablesDict[modifiedName] = "{0}".format(
modifiedName, rawSpecies["compartment"]
)
else:
# ASS - Is this not supposed to be the version without compartments?
self.obs_names.append(modifiedName)
self.obs_map[rawSpecies["identifier"]] = modifiedName
observablesText.append(
"Species {0} {1} #{2}".format(
modifiedName, tmp, rawSpecies["name"]
)
)
observablesDict[modifiedName] = "{0}".format(modifiedName)
speciesTranslationDict[rawSpecies["identifier"]] = tmp
# add the observable in the model
obs_obj = self.bngModel.make_observable()
obs_obj.parse_raw(rawSpecies)
obs_obj.Id = modifiedName
self.bngModel.add_observable(obs_obj)
# TODO: make sure this is replicated in bngModel
sorted(rawSpeciesName, key=len)
for species in rawSpeciesName:
if (
translator[species].getSize() == 1
and translator[species].molecules[0].name not in names
):
names.append(translator[species].molecules[0].name)
mtext = translator[species].str2()
moleculesText.append(mtext)
molec_obj = self.bngModel.make_molecule()
molec_obj.Id = species
# TODO: Make sure we need str2 and not
# just str
if not str(translator[species]) in self.bngModel.molecules:
molec_obj.name = str(translator[species])
else:
molec_obj.name = translator[species].str2()
self.bngModel.add_molecule(molec_obj)
annotationInfo["species"] = speciesAnnotationInfo
self.speciesMemory = []
return (
list(set(moleculesText)),
speciesText,
observablesText,
speciesTranslationDict,
observablesDict,
annotationInfo,
)
[docs] def getInitialAssignments(
self, translator, param, zparam, molecules, initialConditions
):
"""
process the get initial assignments section. This can be used to initialize
parameters or species, so we have to account for both checking both arrays
"""
param2 = param
zparam2 = zparam
initialConditions2 = initialConditions
pparam = {}
initCondMap = {}
for element in param:
pparam[element.split(" ")[0]] = (element.split(" ")[1], None)
for element in zparam:
pparam[element] = (0, None)
for species in self.model.getListOfSpecies():
tmp = self.getRawSpecies(species)
name = tmp["returnID"]
constant = (
"$" if (species.getConstant() or species.getBoundaryCondition()) else ""
)
if name in translator:
if self.noCompartment:
extendedStr = "{1}{0}".format(translator[name], constant)
else:
extendedStr = "@{0}:{2}{1}".format(
species.getCompartment(), translator[name], constant
)
else:
std_name = standardizeName(tmp["name"]) + "()"
# ASS - deal with no compartment case
if self.noCompartment:
extendedStr = "{1}{0}".format(std_name, constant)
else:
extendedStr = "@{0}:{2}{1}".format(
tmp["compartment"], std_name, constant
)
initConc = (
species.getInitialConcentration()
if species.isSetInitialConcentration()
else species.getInitialAmount()
)
pparam[species.getId()] = (initConc, extendedStr)
from copy import copy
for initialAssignment in self.model.getListOfInitialAssignments():
symbol = initialAssignment.getSymbol()
math = libsbml.formulaToString(initialAssignment.getMath())
for element in pparam:
if math is None:
logMess(
"ERROR:SIM210",
f"Initial assignment for '{element}' has no math defined",
)
raise TranslationException(
f"ERROR:SIM210: Initial assignment for '{element}' has no math defined"
)
if element in math:
math = re.sub(
r"(\W|^)({0})(\W|$)".format(element),
r"\1({0})\3".format(pparam[element][0]),
math,
)
# removing non bngl math elements for their equivalents
math = writer.bnglFunction(math, "", []).split(" = ")[1]
param2 = [x for x in param if "{0} ".format(symbol) not in x]
zparam2 = [x for x in zparam if symbol != x]
"""
if (len(param2) != len(param)) or (len(zparam) != len(zparam2)):
param2.append('{0} {1}'.format(symbol,math))
param = param2
zparam = zparam2
"""
try:
if pparam[symbol][1] == None:
param2.append("{0} {1}".format(symbol, math))
param = param2
zparam = zparam2
else:
initialConditions2 = [
x for x in initialConditions if "#{0}".format(symbol) not in x
]
initialConditions2.append(
"{0} {1} #{2}".format(pparam[symbol][1], math, symbol)
)
initCondMap[
"{0} {1} #{2}".format(pparam[symbol][1], math, symbol)
] = symbol
initialConditions = initialConditions2
except:
logMess(
"WARNING:RATE006",
f"Can't find what {symbol} is referring to, some initial conditions might be incorrect.",
)
continue
return param, zparam, initialConditions, initCondMap
[docs] def getSpeciesAnnotation(self):
if self.speciesAnnotation:
return self.speciesAnnotation
self.speciesAnnotation = defaultdict(list)
for species in self.model.getListOfSpecies():
rawSpecies = self.getRawSpecies(species, logEntries=False)
annotationXML = species.getAnnotation()
lista = libsbml.CVTermList()
libsbml.RDFAnnotationParser.parseRDFAnnotation(annotationXML, lista)
if lista.getSize() == 0:
self.speciesAnnotation[rawSpecies["returnID"]] = []
else:
for idx in range(lista.getSize()):
self.speciesAnnotation[rawSpecies["returnID"]].append(
lista.get(idx).getResources()
)
return self.speciesAnnotation
[docs] def getFullAnnotation(self):
if self.speciesAnnotationDict:
return self.speciesAnnotationDict
self.speciesAnnotationDict = defaultdict(lambda: defaultdict(list))
for species in self.model.getListOfSpecies():
rawSpecies = self.getRawSpecies(species, logEntries=False)
annotationXML = species.getAnnotation()
lista = libsbml.CVTermList()
libsbml.RDFAnnotationParser.parseRDFAnnotation(annotationXML, lista)
if lista.getSize() == 0:
continue
else:
for idx in range(lista.getSize()):
for idx2 in range(0, lista.get(idx).getResources().getLength()):
resource = lista.get(idx).getResources().getValue(idx2)
qualifierType = lista.get(idx).getQualifierType()
qualifierDescription = (
bioqual[lista.get(idx).getBiologicalQualifierType()]
if qualifierType
else modqual[lista.get(idx).getModelQualifierType()]
)
self.speciesAnnotationDict[rawSpecies["returnID"]][
qualifierDescription
].append(resource)
return self.speciesAnnotationDict
[docs] def getModelAnnotation(self):
modelAnnotation = []
annotationXML = self.model.getAnnotation()
lista = libsbml.CVTermList()
libsbml.RDFAnnotationParser.parseRDFAnnotation(annotationXML, lista)
if lista.getSize() == 0:
modelAnnotations = []
else:
tempDict = {}
for element in [2, 3, 4, 5, 6]:
if lista.get(element) == None:
continue
tempDict[lista.get(element).getBiologicalQualifierType()] = lista.get(
element
)
if 3 in tempDict:
modelAnnotation = tempDict[3].getResources()
elif 0 in tempDict and (
"GO" in tempDict[0].getResources().getValue(1)
or "kegg" in tempDict[0].getResources().getValue(1)
):
modelAnnotation = tempDict[0].getResources()
elif 5 in tempDict:
modelAnnotation = tempDict[5].getResources()
else:
if lista.get(3) != None and (
"GO" in lista.get(3).getResources().getValue(0)
or "kegg" in lista.get(3).getResources().getValue(0)
):
modelAnnotation = lista.get(3).getResources()
elif lista.get(4) != None and (
"GO" in lista.get(4).getResources().getValue(0)
or "kegg" in lista.get(4).getResources().getValue(0)
):
modelAnnotation = lista.get(4).getResources()
elif lista.get(5) != None and (
"GO" in lista.get(5).getResources().getValue(0)
or "kegg" in lista.get(5).getResources().getValue(0)
):
modelAnnotation = lista.get(5).getResources()
else:
if lista.get(3) != None and (
"reactome" in lista.get(3).getResources().getValue(0)
):
modelAnnotation = lista.get(3).getResources()
elif lista.get(4) != None and (
"reactome" in lista.get(4).getResources().getValue(0)
):
modelAnnotation = lista.get(4).getResources()
elif lista.get(5) != None and (
"reactome" in lista.get(5).getResources().getValue(0)
):
modelAnnotation = lista.get(5).getResources()
elif lista.get(2) != None:
modelAnnotation = lista.get(2).getResources()
return modelAnnotation
[docs] def getSpeciesInfo(self, name):
return self.getRawSpecies(self.model.getSpecies(name))
[docs] def writeLog(self, translator):
rawSpecies = [self.getRawSpecies(x) for x in self.model.getListOfSpecies()]
log["species"].extend([x[0] for x in rawSpecies if x[0] not in translator])
logString = ""
# species stuff
if len(log["species"]) > 0:
logString += "Species we couldn't recognize:\n"
for element in log["species"]:
logString += "\t%s\n" % element
if len(log["reactions"]) > 0:
logString += "Reactions we couldn't infer more about due to \
insufficient information:"
for element in log["reactions"]:
logString += "\t%s + %s -> %s\n" % (
element[0][0],
element[0][1],
element[1],
)
return logString
[docs] def getStandardName(self, name):
if name in self.speciesDictionary:
return self.speciesDictionary[name]
return name
[docs]def standardizeName(name):
"""
Remove stuff not used by bngl
"""
name2 = name
sbml2BnglTranslationDict = {
"^": "",
"'": "",
"*": "m",
" ": "_",
"#": "sh",
":": "_",
"α": "a",
"β": "b",
"γ": "g",
" ": "",
"+": "pl",
"/": "_",
":": "_",
"-": "_",
".": "_",
"?": "unkn",
",": "_",
"(": "",
")": "",
"[": "",
"]": "",
# "(": "__",
# ")": "__",
# "[": "__",
# "]": "__",
">": "_",
"<": "_",
}
for element in sbml2BnglTranslationDict:
name = name.replace(element, sbml2BnglTranslationDict[element])
name = re.sub("[\W]", "", name)
return name