RuleWorld / bionetgen

Rule-based modeling framework
https://bionetgen.org/
MIT License
58 stars 25 forks source link

SBML export won't include compartment information in rate laws #126

Open GoogleCodeExporter opened 9 years ago

GoogleCodeExporter commented 9 years ago
What steps will reproduce the problem?
1. Create a cBNGL model
2. Include the writeSBML() action

What is the expected output? What do you see instead?
Given a reaction @c1:S1 + @c1:S2 -> @c1:S3 k1 , the sbml rate law should be

S1*S2*k1*c1

Instead the compartment is not included in the exported rate law

Original issue reported on code.google.com by jjta...@gmail.com on 2 Jul 2014 at 10:30

GoogleCodeExporter commented 9 years ago
This might require a little more discussion. In bidirectional reactions taken 
place in a given compartment 'c1'

s1 + s2 <-> s3 k1,k2, an sbml model will define the rate function as

((s1*s2*k1) - (s3*k2)) * c1

While in a net file bngl will apply a correction factor of 1/c1 to the forward 
rate, but not to the reverse rate. 

Original comment by jjta...@gmail.com on 2 Jul 2014 at 10:41

GoogleCodeExporter commented 9 years ago
Can you provide a link or a page number in the SBML documentation where they 
discuss this? 

In cBNGL, k1 and k2 have different units, which is why they're treated 
differently. k1 is in units of volume/time, whereas k2 is in units of 1/time. 
That's why k1 is divided by the volume of the compartment (the species 
populations are in units of molecule number). This scaling factor is not 
included in the SBML export, that's true. I'm not sure if this was done on 
purpose or if it's an oversight. Is the 'size' attribute for compartments in 
SBML linked to the kineticLaw at all? If so, this might be why the scaling 
factor is not included in the exported kineticLaw. If not, then we should 
probably fix it so that it is.

Original comment by lh64@cornell.edu on 2 Jul 2014 at 11:03

GoogleCodeExporter commented 9 years ago
If we include the scaling factor in the exported kineticLaw, we may then need 
to change all of the compartment volumes to 1. This may not be a good idea.

Original comment by lh64@cornell.edu on 2 Jul 2014 at 11:06

GoogleCodeExporter commented 9 years ago
Extensive units(counts). cbngl requires those for proper volume correction. 
After double checking the sbml export deals with this by setting the 
hasOnlySubstanceUnits species flag as false (sbml way of indicating molecule 
counts)

>>In cBNGL, k1 and k2 have different units, which is why they're treated 
differently. k1 is in units of volume/time, whereas k2 is in units of 1/time. 
That's why k1 is divided by the volume of the compartment (the species 
populations are in units of molecule number)

Yes, sure. What I actually wanted to say I guess is that SBML actually expects 
a compartment size factor in the rate law for both sides of the equation (I 
discovered this when importing an bng-sbml-exported model in COPASI. So many 
conversions), so some unit conversions might be needed in between during the 
conversion process.

>>I'm not sure if this was done on purpose or if it's an oversight. Is the 
'size' attribute for compartments in SBML linked to the kineticLaw at all? 

No it's not. I checked this manually while importing the models in COPASI. And 
COPASI will complain that your rate law is not true mass action until you 
include the compartment factor.

Original comment by jjta...@gmail.com on 2 Jul 2014 at 11:11

GoogleCodeExporter commented 9 years ago
>>If we include the scaling factor in the exported kineticLaw, we may then need 
to change all of the compartment volumes to 1. This may not be a good idea.

Not at all. You can include the compartment definition on one side, and then 
make a reference to that side on the kinetic law just as when referencing any 
other variable

Original comment by jjta...@gmail.com on 2 Jul 2014 at 11:12

GoogleCodeExporter commented 9 years ago
Here's the Biomodels version of Biomodels 1 ( a fairly small model) for 
illustration purposes. Search for kineticLaw. Here's an example from the first 
reaction

        <kineticLaw metaid="_323346" sboTerm="SBO:0000080">
          <notes>
            <body xmlns="http://www.w3.org/1999/xhtml">
                    <p>kf_0 * B - kr_0 * BL</p>
                        </body>

          </notes>
          <math xmlns="http://www.w3.org/1998/Math/MathML">          
            <apply>
              <times/>
              <ci> comp1 </ci>
              <apply>
                <minus/>
                <apply>
                  <times/>
                  <ci> kf_0 </ci>
                  <ci> B </ci>
                </apply>
                <apply>
                  <times/>
                  <ci> kr_0 </ci>
                  <ci> BL </ci>
                </apply>
              </apply>
            </apply>
          </math>
        </kineticLaw>

Where compartment size is explicitly included on both sides of the equation

Original comment by jjta...@gmail.com on 2 Jul 2014 at 11:24

Attachments:

GoogleCodeExporter commented 9 years ago
Regarding whether the compartment size must be included in an SBML kinetic law, 
the following links expand on this matter

http://sbml.org/Community/Wiki/SBML_Level_3_Core/Reaction_changes/Addition_of_op
tional_compartment_to_reactions

For level2v4, it is referenced in page 72,73, and specially 74 of the 
specification PDF. 

Original comment by jjta...@gmail.com on 3 Jul 2014 at 12:35

GoogleCodeExporter commented 9 years ago
Ok, so it seems to me that what we need to do is set hasOnlySubstanceUnits=true 
when exporting a cBNGL model (i.e., if CompartmentList::getNumCompartments>0) 
and then include the scaling factor in the kineticLaw. It's included in the NET 
file already so it should be easy to do the same thing in writeSBML. This 
follows the "discrete" model example in Sec. 7.3 (page 106) of the L2V4 spec 
(attached).

One thing to consider is that when writing the NET file BNG includes the 
scaling factor as a number since the compartment volume is not stored as a 
separate parameter. However, in the SBML file the 'size' attribute of the 
compartment can be referenced in the kineticLaw simply by including the 
compartment name in the MathML. So I think that's probably what we should do. 
It's simple to do and would mean that the model would be usable in cases where 
the volume changes in time. Obviously, we don't support dynamic volumes right 
now, but hopefully we will at some point in the future.

Does this sound right to you?

Original comment by lh64@cornell.edu on 3 Jul 2014 at 6:26

Attachments:

GoogleCodeExporter commented 9 years ago
Yes, sounds perfect.

Original comment by jjta...@gmail.com on 3 Jul 2014 at 4:35