BioNetGen does not export "initialAssigments" in SBML files

NikosMemmos commented 1 year ago

Hi everyone,

I have been using BioNetGen to export my model in an SBML file to perform Global Sensitivity Analysis, using GloabalSensitivityAnalysis package in Julia. One issue I have been encountering is that some of the parameters in my models are the initial conditions for certain species, but I observed that changing the parameter, does not change the initial condition automatically during the sensitivity analysis (see topic here: https://discourse.julialang.org/t/performing-global-sensitivity-analysis-on-a-big-reaction-network-imported-as-a-sbml-file/93000)

After talking with @paulflang who has been developing the SBMLToolkit, he observed that BioNetGen does not export initialAssigments but plugs in the numeric values in the SBML file, which causes this issue.

Is there any way to fix the problem? Thanks for your time.

paulflang commented 1 year ago

So IIUC what happens is that if you do A_0 10.0 in the parameter section and then A() A_0 when you seed the species, then the parameter A_0 does not appear in the exported SBML. However, it would be good if it appeared in the initialAssignments, such that things like parameter estimation of A_0 can be done.

ASinanSaglam commented 1 year ago

One potential work around to this problem is to write expressions that define starting values for species which are then exported as initialAssignments that can be manipulated for downstream parameter estimation. However, there is a bug here and you have to a small change in the SBML itself

Here is an example BNGL

begin parameters
  a_cnt 100
end parameters

begin species
 A(b) a_cnt
 B(a) a_cnt*10
end species

in SBML this gets exported to (removing some extra stuff)

    <listOfSpecies>
      <species id="S1" compartment="cell" initialAmount="100" name="A(b)"/>
      <species id="S2" compartment="cell" initialAmount="1000" name="B(a)"/>
    </listOfSpecies>
    <listOfParameters>
      <!-- Independent variables -->
      <parameter id="a_cnt" value="100"/>
      <!-- Dependent variables -->
      <parameter id="_InitialConc1" constant="true"/>
    </listOfParameters>
    <listOfInitialAssignments>
      <!-- Dependent variables -->
      <initialAssignment symbol="_InitialConc1">
        <math xmlns="http://www.w3.org/1998/Math/MathML">
          <apply>
            <times/>
            <ci> a_cnt </ci>
            <cn> 10 </cn>
          </apply>
        </math>
      </initialAssignment>
    </listOfInitialAssignments>

if you just remove the intermediary parameter that's automatically generated and replace the symbol the initial assignment is adjusting to the correct species, this will work as expected. E.g.

    <listOfSpecies>
      <species id="S1" compartment="cell" initialAmount="100" name="A(b)"/>
      <species id="S2" compartment="cell" initialAmount="1000" name="B(a)"/>
    </listOfSpecies>
    <listOfParameters>
      <!-- Independent variables -->
      <parameter id="a_cnt" value="100"/>
    </listOfParameters>
    <listOfInitialAssignments>
      <!-- Dependent variables -->
      <initialAssignment symbol="S2">
        <math xmlns="http://www.w3.org/1998/Math/MathML">
          <apply>
            <times/>
            <ci> a_cnt </ci>
            <cn> 10 </cn>
          </apply>
        </math>
      </initialAssignment>
    </listOfInitialAssignments>

I'll try to take a look at the code today and see if I can properly automate this.

paulflang commented 1 year ago

Thanks a lot for looking into this, @ASinanSaglam . It would be indeed good to automate this, cause users would hardly know what to do and even if they did, for large models manual changes can be very cumbersome and error prone. However, IIUC you also need to add S1 to the initial assignments, probably roughly like so (otherwise you can change a_cnt, without anything happening to the initial condition of S1 in SBML).

    <listOfInitialAssignments>
      <!-- Dependent variables -->
      <initialAssignment symbol="S1">
        <math xmlns="http://www.w3.org/1998/Math/MathML">
            <ci> a_cnt </ci>
        </math>
      </initialAssignment>
      <initialAssignment symbol="S2">
        <math xmlns="http://www.w3.org/1998/Math/MathML">
          <apply>
            <times/>
            <ci> a_cnt </ci>
            <cn> 10 </cn>
          </apply>
        </math>
      </initialAssignment>
    </listOfInitialAssignments>

ASinanSaglam commented 1 year ago

I think I got something working. Before I merge, I wanted to ask, @paulflang does this SBML block look about right to you?

    <listOfInitialAssignments>
      <!-- All initial values -->
      <initialAssignment symbol="S1">
        <math xmlns="http://www.w3.org/1998/Math/MathML">
          <ci> a_init </ci>
        </math>
      </initialAssignment>
      <initialAssignment symbol="S2">
        <math xmlns="http://www.w3.org/1998/Math/MathML">
          <apply>
            <times/>
            <cn> 15 </cn>
            <ci> k1 </ci>
          </apply>
        </math>
      </initialAssignment>
      <initialAssignment symbol="S3">
        <math xmlns="http://www.w3.org/1998/Math/MathML">
          <ci> p_test </ci>
        </math>
      </initialAssignment>
      <initialAssignment symbol="S4">
        <math xmlns="http://www.w3.org/1998/Math/MathML">
          <ci> 10 </ci>
        </math>
      </initialAssignment>
      <!-- Dependent variables -->
      <initialAssignment symbol="p_test">
        <math xmlns="http://www.w3.org/1998/Math/MathML">
          <apply>
            <times/>
            <apply>
              <times/>
              <cn> 10 </cn>
              <ci> k1 </ci>
            </apply>
            <ci> factor </ci>
          </apply>
        </math>
      </initialAssignment>
    </listOfInitialAssignments>

input BNGL looks something like this

begin model

begin parameters
  a_init 10
  k1 1
  k2 1
  factor 5
  p_test 10*k1*factor
end parameters

begin molecule types
 A(b)
 B(a)
 C()
end molecule types

begin observables
 Molecules A_free A(b)
end observables

begin species
 A(b) a_init
 A(b!1).B(a!1) 15*k1
 B(a) p_test
 C() 10
end species

begin reaction rules
 A(b) + B(a) <-> A(b!1).B(a!1) k1,k2
end reaction rules

end model

Let me know if I'm missing something, otherwise I'll merge this and it'll be available with the next version of BNG. The new release will be done in a couple weeks.

paulflang commented 1 year ago

Thanks @ASinanSaglam , looks good to me, except of

      <initialAssignment symbol="S4">
        <math xmlns="http://www.w3.org/1998/Math/MathML">
          <ci> 10 </ci>
        </math>
      </initialAssignment>

which is not needed (the initialConcentration or initialAmount field of species S4, should simply be set to 10. No need to override it with an initialAssignment).

But please, try to import the exported SBML with a tool of your choice (e.g. Copasi) to check for any warnings/errors, before merging the PR. Some things are hard to pick up by eye.

ASinanSaglam commented 1 year ago

I checked it via COPASI, I just missed that an integer value should be a type cn and not ci in MathML so that's fixed. Now my basic test gets imported into COPASI without any issues.

Is there a major issue for exporting all species in initial assignments that I'm missing? I know exporting scalar values in initial assignments block is redundant and will make it less human readable but is there a technical issue with it?

I'd like to keep it consistent and export all species in initial assignments, that ensures that all initial values are stored in a single location instead of having them spread in multiple places and just in case that's something some other library likes to have.

paulflang commented 12 months ago

Sorry, only see this now.

Is there a major issue for exporting all species in initial assignments that I'm missing?

It is correct, but I think it is

less human readable, as you say
an unconventional thing to do
several current or future software packages may not support SBML with initial assignments, especially if the developers focussed on things other than compliance with a lot of different nuances of SBML. Given the importance of BNGL in the systems biology software ecosystem, it would be a shame if models that would have previously been exported/imported easily would no longer be exported/imported due to this bugfix.

ensures that all initial values are stored in a single location instead of having them spread in multiple places

I think the idea is that initial values are defined with the species per default, but can be overridden in initialAssigments. So there are essentially just two places to look at, the species definition and the (often empty) list of initialAssignments.

So if it is an Int or Float, I would not make it an initialAssignment. But I think what you are doing is technically correct, and it is ultimately your choice how to implement it.

ASinanSaglam commented 12 months ago

No worries at all, I just wanted to push a version of the fix, it is not released yet. We should have a new release early next week.

While I'm not the biggest fan of splitting where the values are stored, I agree it'll make it more confusing and be potentially backwards incompatible. I will change it so that the numerical values are skipped for the initialAssignments block and adhere to SBML conventions.

RuleWorld / bionetgen

BioNetGen does not export "initialAssigments" in SBML files #264