Exporting to other formats

[himoto]

Hi @formersbach, I would like to discuss how to export biomass models to other formats.
Here you can share ideas and code to achieve this.
Thank you!

[formersbach]

[formersbach]

That's an example of how the ODE file for a simple E + S -> ES -> P system would look like. Accessing for example dsdt would return just the equation. Accessing dydt before substitution would return the entire dynamic system. Returning list(dydt) is equivalent to the returned value as it is currently implemented in Biomass, with minor differences due to sympy calculating less digits per default.
I have some questions/problems regarding the actual implementation, which I will formulate in separate comments!

[himoto]

Last two lines can be replaced with:

_dydt = sympy.Matrix(...)
all_subs = {}
all_subs.update({sympy.Symbol(pname): eval(f"x[C.{pname}]") for pname in C.NAMES})
all_subs.update({sympy.Symbol(sname): eval(f"y[V.{sname}]") for sname in V.NAMES})
dydt = _dydt.subs(all_subs)

[himoto]

Correction:

_dydt = sympy.Matrix(...)
all_subs = {}
all_subs.update({sympy.Symbol(pname): eval(f"{x[i]}") for i, pname in enumerate(C.NAMES)})
all_subs.update({sympy.Symbol(sname): eval(f"y[V.{sname}]") for sname in V.NAMES})
dydt = _dydt.subs(all_subs)

[formersbach]

As you can see in the image each symbol that is used in the expressions needs to be individually declared. This is impractical for large models with > 10 params/species. The only approach I've come up with so far is modifying globals, which as we've discussed before is not an optimal solution.
Could the species/parameter dataclass be modified to make that possible? Any ideas?

[formersbach]

Can I pass the flux equations as symbolic equations to the ODE construction or are the flux equations evaluated outside of ode.py?

[himoto]

I created the kinetics branch, in which users can access information about kinetics.
Suppose that all models are generated from Text2Model.

1. Prepare a text file describing Michaelis-Menten kinetics (michaelis_menten.txt)

E + S <--> ES | kf=0.003, kr=0.001 | E=100, S=50
ES --> E + P | kf=0.002

@obs Substrate: u[S]
@obs E_free: u[E]
@obs E_total: u[E] + u[ES]
@obs Product: u[P]
@obs Complex: u[ES]

@sim tspan: [0, 100]

2. Convert text into an executable BioMASS model

from biomass import Text2Model
description = Text2Model("michaelis_menten.txt")
description.convert()

Now we can access information about kinetics via description.kinetics:

for r in description.kinetics:
    print("reactants:", r.reactants)
    print("products:", r.products)
    print("rate:", r.rate)
    print("--------")

reactants: ('E', 'S')
products: ('ES',)
rate: kf1 * E * S
--------
reactants: ('ES',)
products: ('E', 'S')
rate: kr1 * ES
--------
reactants: ('ES',)
products: ('E', 'P')
rate: kf2 * ES
--------

Do you think this is useful for exporting models to other formats?

[formersbach]

Wow that's super useful!! Thanks a lot!
There are two caveats I noticed:

1. It's only available for the Text2Model object, which means the model could only be exported to sbml during creation and not afterwards.
2. The rates in .kinetics do not match up with the rates in reaction_network.py. For example in the menten model there are three fluxes/rates while in reaction_network.py only two rates are present. As far as I checked bidirectional reactions are treated differently for self.reaction and self.kinetics. Is that intended?

[himoto]

1. Yes, currently this is only available for the Text2Model object. If you think the kinetics info is useful, I am planning to copy it in biomass.ModelObject class as well. Please note that we always require models to be created from Text2Model.
2. Yes, the resulting kinetics do not match up with the flux vector in reaction_network.py as I made each reaction unidirectional. I think for me, this is ok for exporting.

[formersbach]

Concerning 1.
Awesome that would be really nice!
Concerning 2.
The problem with keeping each reaction unidirectional is that the fluxes are not synced with the ODE's anymore. Therefore I cannot generate the full ODE correctly, since the fluxes cannot be properly substituted. In the menten example:

ODE:
dydt[V.E] = - v[1] + v[2]
Flux:
v[1] = x[C.kf1] * y[V.E] * y[V.S] - x[C.kr1] * y[V.ES]
v[2] = x[C.kf2] * y[V.ES]
Kinetics:
k[1] = kf1 * E * S
k[2] = kr1 * ES
k[3] = kf2 * ES

If I now try to substitute the ODE I'll get dedt = - kf1 * E * S + kr1 * ES.
Hope the example clarifies the problem ^^

[himoto]

Since we can easily access the information about reactants and products in each reaction, we can create the model's stoichiometry matrix (S). Please refer to the latest updates in the kinetics branch.

>>> description.stoichiometry_matrix.toarray()
array([[-1,  1,  1],
       [-1,  1,  0],
       [ 1, -1, -1],
       [ 0,  0,  1]])

dydt = S * v, where S is stoichiometry matrix and v is reaction vector.
Therefore,
d[E]/dt = -v1 + v2 + v3 = - kf1 * E * S + kr1 * ES + kf2 * ES
d[S]/dt = -v1 + v2 = - kf1 * E * S + kr1 * ES
d[ES]/dt = v1 - v2 - v3 = kf1 * E * S - kr1 * ES - kf2 * ES
d[P]/dt = v3 = kf2 * ES

[formersbach]

Great! That's even more useful than having the reaction rates in the same format!
Now we have all the building blocks for the entire system right?

[formersbach]

With the changes you implemented in the kinetics branch it's relatively straightforward to turn models into SBML now. Thanks a lot! However there are further decisions that we need to make about the implementation:

1. Where should the sbml conversion code be added? Imo adding it to the model object would be the most sensible option. It's approximately 100 lines of code. Being able to call model.2sbml(options) and model.safe_sbml(options) is where I would expect it as a user.
2. What parts of the model should be transferred into SBML? As you likely are aware SBML supports a sheer infinite amount of customization. Currently the SBML file contains: Species, Reactions, Initial values, parameter values. Currently the things that could be implemented (off the top of my head) are:
   • parameter ranges
   • unit specifications
   • different conditions
   • events
   • observables
   • parameter optimization conditions
   • ....

Please let me know what you think should be in the sbml file and what not! You probably also know customization options in Biomass that I don't so I'd really appreciate your input.

[himoto]

1. Please create export.py in biomass/. I think we should start from exporting model from Text2Model object. Therefore, if you create ExportModel class in the module, if should receive Text2Model class.

# Example
class ExportModel(object):

def __init__(self, description: Text2Model):
    self.description = description

to_sbml(self, ...):
    ...

2. I think if would be nice to start from exporting minimal information, which we can extract from Text2Model class. If we find it works, we can move on to the next step, i.e., including more information.