6 # -*- coding: utf-8 -*-
"""
Created on Mon Jun 17 11:19:37 2013
@author: proto
"""
import libsbml
import json
from optparse import OptionParser
[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))
[docs]class SBML2JSON:
def __init__(self, model):
self.model = model
self.getUnits()
self.moleculeData = {}
[docs] def getUnits(self):
self.unitDictionary = {}
# unitDictionary['substance'] = [libsbml.UNIT_KIND_MOLE,1]
# unitDictionary['volume'] = [libsbml.UNIT_KIND_LITER,1]
# unitDictionary['area'] = [libsbml.UNIT_KIND_METER,1]
# unitDictionary['length'] = [libsbml.UNIT_KIND_METER,1]
# unitDictionary['time'] = [libsbml.UNIT_KIND_SECOND,1]
for unitDefinition in self.model.getListOfUnitDefinitions():
unitList = []
for unit in unitDefinition.getListOfUnits():
unitList.append([unit.getKind(), unit.getScale(), unit.getExponent()])
self.unitDictionary[unitDefinition.getId()] = unitList
"""
Table 3: SBML's built-in units.
Name Possible Scalable Units Default Units
substance mole, item mole
volume litre, cubic metre litre
area square metre square metre
length metre metre
time second second
"""
[docs] def getParameters(self):
parameters = []
prx = {
"name": "Nav",
"value": "6.022e8",
"unit": "",
"type": "Avogadro number for 1 um^3",
}
parameters.append(prx)
for parameter in self.model.getListOfParameters():
parameterSpecs = {
"name": parameter.getId(),
"value": parameter.getValue(),
"unit": parameter.getUnits(),
"type": "",
}
"""
if parameterSpecs[0] == 'e':
parameterSpecs = ('are',parameterSpecs[1])
if parameterSpecs[1] == 0:
zparam.append(parameterSpecs[0])
"""
if parameter.getUnits() in self.unitDictionary:
for factor in self.unitDictionary[parameter.getUnits()]:
parameterSpecs["value"] *= 10 ** (factor[1] * factor[2])
parameterSpecs["unit"] = "{0}*1e{1}".format(
parameterSpecs["unit"], factor[1] * factor[2]
)
if (
"mole" in parameter.getUnits()
and "per_mole" not in parameter.getUnits()
):
parameterSpecs["value"] *= float(6.022e8)
parameterSpecs["unit"] = "{0}*{1}".format(
parameterSpecs["unit"], "avo.num"
)
# if parameter.getUnits() == '':
# parameterSpecs['value'] *= float(6.022e8*1000)
# parameterSpecs['unit'] = '{0}*{1}'.format(parameterSpecs['unit'],'avo.num*1000')
parameters.append(parameterSpecs)
prx = {"name": "rxn_layer_t", "value": "0.01", "unit": "um", "type": ""}
ph = {"name": "h", "value": "rxn_layer_t", "unit": "um", "type": ""}
pRs = {"name": "Rs", "value": "0.002564", "unit": "um", "type": ""}
pRc = {"name": "Rc", "value": "0.0015", "unit": "um", "type": ""}
parameters.append(prx)
parameters.append(ph)
parameters.append(pRs)
parameters.append(pRc)
parameterDict = {idx + 1: x for idx, x in enumerate(parameters)}
return parameterDict
def __getRawCompartments(self):
"""
extracts information about the compartments in a model
*returns* name,dimensions,size
"""
compartmentList = {}
for compartment in self.model.getListOfCompartments():
name = compartment.getId()
size = compartment.getSize()
outside = compartment.getOutside()
dimensions = compartment.getSpatialDimensions()
compartmentList[name] = [dimensions, size, outside]
return compartmentList
[docs] def getOutsideInsideCompartment(self, compartmentList, compartment):
"""
Gets the containing compartment for this compartment
"""
outside = compartmentList[compartment][2]
for comp in compartmentList:
if compartmentList[comp][2] == compartment:
return outside, comp
return outside, -1
[docs] def getMolecules(self):
"""
*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
"""
compartmentList = self.__getRawCompartments()
molecules = []
release = []
for idx, species in enumerate(self.model.getListOfSpecies()):
compartment = species.getCompartment()
if compartmentList[compartment][0] == 3:
typeD = "3D"
outside, inside = self.getOutsideInsideCompartment(
compartmentList, compartment
)
diffusion = "KB*T/(6*PI*mu_{0}*Rs)".format(compartment)
else:
typeD = "2D"
diffusion = "KB*T*LOG((mu_{0}*h/(SQRT(4)*Rc*(mu_{1}+mu_{2})/2))-gamma)/(4*PI*mu_{0}*h)".format(
compartment, outside, inside
)
self.moleculeData[species.getId()] = [compartmentList[compartment][0]]
moleculeSpecs = {
"name": species.getId(),
"type": typeD,
"extendedName": species.getName(),
"dif": diffusion,
}
initialConcentration = species.getInitialConcentration()
if initialConcentration == 0:
initialConcentration = species.getInitialAmount()
if species.getSubstanceUnits() in self.unitDictionary:
for factor in self.unitDictionary[species.getSubstanceUnits()]:
initialConcentration *= 10 ** (factor[1] * factor[2])
if "mole" in species.getSubstanceUnits():
initialConcentration /= float(6.022e8)
if species.getSubstanceUnits() == "":
initialConcentration /= float(6.022e8)
isConstant = species.getConstant()
# isBoundary = species.getBoundaryCondition()
if initialConcentration != 0:
if compartmentList[compartment][0] == 2:
objectExpr = "{0}[{1}]".format(inside.upper(), compartment.upper())
else:
objectExpr = "{0}".format(compartment)
releaseSpecs = {
"name": "Release_Site_s{0}".format(idx + 1),
"molecule": species.getId(),
"shape": "OBJECT",
"quantity_type": "NUMBER_TO_RELEASE",
"quantity_expr": initialConcentration,
"object_expr": objectExpr,
}
release.append(releaseSpecs)
# self.speciesDictionary[identifier] = standardizeName(name)
# returnID = identifier if self.useID else \
molecules.append(moleculeSpecs)
# self.sp eciesDictionary[identifier]
moleculesDict = {idx + 1: x for idx, x in enumerate(molecules)}
releaseDict = {idx + 1: x for idx, x in enumerate(release)}
return moleculesDict, releaseDict
[docs] def getPrunnedTree(self, math, remainderPatterns):
"""
remove mass action factors
"""
while (math.getCharacter() == "*" or math.getCharacter() == "/") and len(
remainderPatterns
) > 0:
if libsbml.formulaToString(math.getLeftChild()) in remainderPatterns:
remainderPatterns.remove(libsbml.formulaToString(math.getLeftChild()))
math = math.getRightChild()
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 getInstanceRate(self, math, compartmentList, reversible, rReactant, rProduct):
# remove compartments from expression
math = self.getPrunnedTree(math, compartmentList)
if reversible:
if math.getCharacter() == "-" and math.getNumChildren() > 1:
rateL, nl = self.removeFactorFromMath(
math.getLeftChild().deepCopy(), rReactant, rProduct
)
rateR, nr = self.removeFactorFromMath(
math.getRightChild().deepCopy(), rProduct, rReactant
)
else:
rateL, nl = self.removeFactorFromMath(math, rReactant, rProduct)
rateL = "if({0} >= 0 ,{0},0)".format(rateL)
rateR, nr = self.removeFactorFromMath(math, rReactant, rProduct)
rateR = "if({0} < 0 ,-({0}),0)".format(rateR)
nl, nr = 1, 1
else:
rateL, nl = self.removeFactorFromMath(math.deepCopy(), rReactant, rProduct)
rateR, nr = "0", "-1"
if reversible:
pass
return rateL, rateR
[docs] def removeFactorFromMath(self, math, reactants, products):
remainderPatterns = []
highStoichoiMetryFactor = 1
for x in reactants:
highStoichoiMetryFactor *= factorial(x[1])
y = [i[1] for i in products if i[0] == x[0]]
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)
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]
math = self.getPrunnedTree(math, remainderPatterns)
rateR = libsbml.formulaToString(math)
for element in remainderPatterns:
rateR = "if({0} >0,({1})/{0} ,0)".format(element, rateR)
if highStoichoiMetryFactor != 1:
rateR = "{0}*{1}".format(rateR, int(highStoichoiMetryFactor))
return rateR, math.getNumChildren()
[docs] def adjustParameters(self, stoichoimetry, rate, parameters):
"""
adds avogadros number and other adjusting factors to the reaction rates
"""
for parameter in parameters:
if (
parameters[parameter]["name"] in rate
and parameters[parameter]["unit"] == ""
):
print(parameters[parameter])
if stoichoimetry == 2:
parameters[parameter]["value"] *= float(6.022e8)
parameters[parameter]["unit"] = "Bimolecular * NaV"
elif stoichoimetry == 0:
parameters[parameter]["value"] /= float(6.022e8)
parameters[parameter]["unit"] = "0-order / NaV"
elif stoichoimetry == 1:
parameters[parameter]["unit"] = "Unimolecular"
print(parameters[parameter])
[docs] def getReactions(self, sparameters):
"""
returns a list with reactant,product and fwdRate
"""
reactionSpecs = []
for index, reaction in enumerate(self.model.getListOfReactions()):
reactant = [
(reactant.getSpecies(), reactant.getStoichiometry())
for reactant in reaction.getListOfReactants()
if reactant.getSpecies() != "EmptySet"
]
product = [
(product.getSpecies(), product.getStoichiometry())
for product in reaction.getListOfProducts()
if product.getSpecies() != "EmptySet"
]
kineticLaw = reaction.getKineticLaw()
rReactant = [
(x.getSpecies(), x.getStoichiometry())
for x in reaction.getListOfReactants()
if x.getSpecies() != "EmptySet"
]
rProduct = [
(x.getSpecies(), x.getStoichiometry())
for x in reaction.getListOfProducts()
if x.getSpecies() != "EmptySet"
]
parameters = [
(parameter.getId(), parameter.getValue())
for parameter in kineticLaw.getListOfParameters()
]
math = kineticLaw.getMath()
reversible = reaction.getReversible()
compartmentList = []
for compartment in self.model.getListOfCompartments():
compartmentList.append(compartment.getId())
rateL, rateR = self.getInstanceRate(
math, compartmentList, reversible, rReactant, rProduct
)
# finalReactant = [x[0]]
# testing whether we have volume-surface interactions
rcList = []
prdList = []
for element in reactant:
orientation = (
","
if len(set(self.moleculeData[x[0]][0] for x in reactant)) > 1
and self.moleculeData[element[0]] == "3"
else "'"
)
rcList.append("{0}{1}".format(element[0], orientation))
for element in product:
orientation = (
","
if len(set(self.moleculeData[x[0]][0] for x in reactant)) > 1
and self.moleculeData[element[0]] == "3"
else "'"
)
prdList.append("{0}{1}".format(element[0], orientation))
if rateL != 0:
tmp = {}
tmp["reactants"] = " + ".join(rcList)
tmp["products"] = " + ".join(prdList)
tmp["fwd_rate"] = rateL
reactionSpecs.append(tmp)
if rateR != 0:
tmp = {}
tmp["reactants"] = " + ".join(prdList)
tmp["products"] = " + ".join(rcList)
tmp["fwd_rate"] = rateR
reactionSpecs.append(tmp)
self.adjustParameters(len(reactant), rateL, sparameters)
self.adjustParameters(len(product), rateR, sparameters)
reactionDict = {idx + 1: x for idx, x in enumerate(reactionSpecs)}
return reactionDict
# SBML USE INSTANCE RATE
# HOW TO GET THE DIFFUSION CONSTANT
[docs]def main():
# command line arguments
parser = OptionParser()
parser.add_option(
"-i",
"--input",
dest="input",
default="bngl2mcell/rec_dim_sbml.xml",
type="string",
help="The input SBML file in xml format. Default = 'input.xml'",
metavar="FILE",
)
parser.add_option(
"-o",
"--output",
dest="output",
type="string",
help="the output JSON file. Default = <input>.py",
metavar="FILE",
)
(options, args) = parser.parse_args()
reader = libsbml.SBMLReader()
nameStr = options.input
if options.output == None:
outputFile = nameStr + ".py"
else:
outputFile = options.output
# libsbml initialization stuff
document = reader.readSBMLFromFile(nameStr)
if document.getModel() == None:
print("No such input file")
return
# get data
parser = SBML2JSON(document.getModel())
parameters = parser.getParameters()
molecules, release = parser.getMolecules()
reactions = parser.getReactions(parameters)
definition = {}
definition["par_list"] = parameters
definition["mol_list"] = molecules
definition["rxn_list"] = reactions
definition["rel_list"] = release
print("Writing output to {0}".format(outputFile))
# output
with open(outputFile, "w") as f:
json.dump(definition, f, sort_keys=True, indent=1, separators=(",", ": "))
if __name__ == "__main__":
main()