Using Systems Biology Markup Langague (SBML) in Daphne for Model Sharing

What is SBML?

The systems biology markup language [(SBML)](http://sbml.org/Main_Page) is a popular and widely accepted XML-based, open-source language for representing and exchanging models across analysis and simulation tools. By utilizing the capabilities provided by SBML though [libSBML](http://sbml.org/Software/libSBML), Daphne is able to translate user-defined immunological models into SBML formatted files, which can then be shared, evaluated and cooperatively developed.

Daphne model sharing

The fundamental goal of the MSI model sharing plan is to enable our software platform to allow users to share, evaluate and cooperatively develop models underlying agent-based immune response simulations. For this purpose, we decided to use an information representation standard developed by the systems biology community known as the Systems Biology Markup Language (SBML) [Hucka2003]. The latest generation of SBML (Level 3) is modular, which means that there is a core set of elements for model description and optional packages that layer additional features on top of the core. By utilizing the capabilities provided by SBML, we will enable the representation and exchange of user-defined models.

Our model sharing efforts took on a three-pronged approach: implementation of SBML Level 3 core, implementation of Level 3 packages, and development of new SBML capabilities.

Implementation of SBML Level 3 core models

The first step towards our model sharing goals was to provide the capability to export immunological models as SBML Level 3 (core) models and to be able to import Level 3 (core) models for subsequent model analysis and simulations. As we have now enabled these features in our tool, user-defined models can be translated into SBML files, and can be easily read back into our software platform. Since the SBML Level 3 core package cannot fully encode the models in Daphne, we use Level 3 packages to further specify model information. See next section for more information. When there is no suitable Level 3 package, we use annotation, where possible, to specify model information.

Implementation of SBML Level 3 packages

The second step was to evaluate and, where possible, utilize SBML packages that are available or under development that might extend the Daphne model features that can be encoded without annotation (without flattening). The latest generation of SBML (Level 3) is designed around the concept of modularity, which specifies a core set of data constructs for model encoding as well as optional packages that layer additional features on top of the core. To achieve our modeling goals, we have used a combination of constructs provided by different SBML packages to describe the modeling features used in our immunological simulations. These are the aspects where our model sharing goals have benefited from the recent usage of SBML packages:

  1. Our simulation system is inherently multi-scaled and the presence of autonomous entities (molecules or cells) and their local behaviors (actions and interactions) determine the evolution of the overall system (germinal center). The comp SBML extension, which is already available, is being used to provide the capabilities needed for this agent-based type of simulation where we aggregate modeling components into cell submodels that help encode the topological relationships within a model.
  2. All of the modeled intracellular dynamics and physiological environments of biochemical reactants are spatially distributed. The capability to describe spatial systems is currently available via the SBML spatial package. However, as this package has only been partially implemented and is still under development, we have purely used it to allow users to export spatial models. Importing models to ensure reproducibility will only be achievable once the standard is fully developed and available to the community.
  3. Though structural constructs are currently fixed in SBML, an accurate depiction of immune response in the germinal center requires modeling cellular events such as cell proliferation, differentiation and cell death. In order to deliver this capability, we have preliminarily implemented the constructs in the first draft of the package specification for the proposed dyn extension. See next section for more information.
  4. Groups of reactions in our simulator can be aggregated into reaction complexes, which users can create to refer to specific groups of reactions that as a whole represent a certain process inside the cell. To properly encode this in SBML, we have provided support for the SBML groups package (in development), which allows us to create a collection of reactions and indicate their relation.

Development of SBML capabilities

Our model sharing goals benefited from already available packages, but also required further development of SBML capabilities. For this reason, in the third step we actively participated in the development of new SBML capabilities

Dynamic Structure (dyn) Package

We collaborated with Dr. Chris Myers (University of Utah), in the development of the dyn package, which allows for the description of biological processes such as cell birth, differentiation and cell death in our models. Towards this end, we drafted a proposal that followed the basic architectural principles required (“Development Process for SBML Level 3 - SBML.org,” n.d.). Our proposal reflected the purpose, general approach and orientation of a package which will to allow for the representation of dynamic cellular behavior in multicellular systems. Having obtained favorable feedback after submitting the proposal for community approval, we proceeded to the next stage of development. A copy of the proposal can be obtained at the SBML Dynamic Structures website. http://sbml.org/Documents/Specifications/SBML_Level_3/Packages/Dynamic_Structures_(dyn)

This work was presented at the 2014 COMBINE meeting (Harold Gomez, Bioinformatics Graduate Student in the Kepler lab), a yearly workshop-style event about standards and standardization efforts in computational models, where it received constructive feedback from the SBML community.

Contrary to the proposal stage, the SBML package specification development required significant effort and extensive software implementation. We established a Package Working Group (PWG), a community standard whose goal is to formally bring together individuals who have relevant expertise to collaborate throughout the life of the package. We have also began to specify new encoding constructs which are consistent with the dyn package proposal, and plan to provide software implementations for these newly defined constructs in two software platforms in the near future.

With our group taking the lead, a proposal that follows the basic architectural principles required (“Development Process for SBML Level 3 - SBML.org,” n.d.) was drafted and accepted by the community. The proposal is available at the following URL:

http://sbml.org/Community/Wiki/SBML_Level_3_Proposals/Dynamic_Structures

Next, we proceeded to draft a package specification. This involved translating the purpose, general approach and orientation of the package- which is to allow for the representation of dynamic cellular behavior in multicellular systems, into construct definitions that people can use in SBML models. It is important to note that key insights about both the development of package specifications in SBML and the proposed constructs of the dyn extension were obtained in the bi-yearly SBML conferences COMBINE and HARMONY. The first draft of the dyn package specification can be found here:

http://sbml.org/Documents/Specifications/SBML_Level_3/Packages/Dynamic_Structures_%28dyn%29

Package implementation in libSBML and JSBML

SBML comes with two free and open-source application programming interfaces (APIs): libSBML and JSBML [Bornstein2008] which provide the programmatic capabilities for reading, writing, manipulating, translating, and validating models expressed in SBML format. Prior to final approval of our extension from SBML editors, there must be at least two software systems that offer manipulation of the dyn constructs. To achieve this, we have translated our draft specification into libSBML and JSBML APIs that serve to implement dyn SBML functionality via programming languages. This draft implementation is already available in the latest release of both libSBML and JSBML.

The next step in the development will be for two tools to implement the dyn package to show the community the usefulness of proposed constructs and to encourage other adopters. This process generally requires a few more iterations of the specification before formal acceptance as a full-fledged SBML package, but people can already start using it as we have provided library implementations.

Development of MOCCASIN

We were part of the core group of developers for MOCCASIN [Rodriguez2015] which stands for “Model ODE Converter for Creating Awesome SBML INteroperability”. MOCCASIN is a Python tool primarily designed to interpret certain basic forms of ODE simulation models written in MATLAB or Octave and translate them into SBML format. The resulting software addresses the needs of a growing body of researchers who find it discouragingly difficult to export their MATLAB/Octave models to open and widely-used formats in systems biology such as SBML.

At its core, MOCCASIN uses an algorithm developed by Fages and colleagues [Fages2014] to infer reaction systems from sets of ODEs. To parse MATLAB files and produce input to the reaction-inference algorithm, MOCCASIN uses a custom MATLAB parser written using PyParsing [McGuire2015], and a variety of post-processing operations to interpret MATLAB file contents; MATLAB isn’t required to run MOCCASIN. Users who want to use MOCCASIN can do it either via a command line or a graphical user interface. For more information on MOCCASIN, please refer to the repository at: https://github.com/sbmlteam/moccasin.

Currently, MOCCASIN is limited to MATLAB inputs in which a model is contained in a single file. Additionally, the file must set up a system of differential equations as a function defined in the file, and make a call to one of the MATLAB odeNN family of solvers (e.g., ode45, ode15s, etc.).

How to use:

To encode models into SBML in Daphne, click on “File” -> “Export in SBML”. Two files will be added to the chosen output directory: an XML file corresponding to the SBML model, and a corresponding log file.

Similarly, to load an SBML model into Daphne for simulation, click on “File” -> “Import from SBML”. Prior to loading, models in selected SBML files are converted to the latest iteration of SBML and a log file is deposited in the source directory.

Current limitations:

The latest generation of SBML, which is Level 3 Version 1, is modular as it has a defined core set of features and [optional packages](http://sbml.org/Community/Wiki) that add additional elements. As an initial step, we’ve only supported core features to describe our models. Model components that could not be encoded in this way were added as tool-specific annotations, which are not exchangeable across tools. Though initially all encoded models were treated as single-celled models and reduced to their underlying set of reactions, this has been expanded to support complex multicellular simulation setups. Dynamic cellular events such cell movement, differentiation, proliferation and cell death are currently excluded.

Future directions:

Our goal is to incrementally support SBML packages to enrich model description as we move forward. Since packages have different strengths and there is no single package that answers all our needs, we will use a combination of constructs provided by different packages (comp, spatial, dyn) to accurately describe the complexities of immunological simulations in Daphne.

Example:

<?xml version='1.0' encoding='UTF-8' standalone='no'?>
<sbml xmlns="http://www.sbml.org/sbml/level3/version1/core" level="3" version="1">
  <model id="ReceptorDynamics" areaUnits="microMetreSquared"
  name="PlazaSur Receptor Dynamics" metaid="_ReceptorDynamics"
  timeUnits="minute" lengthUnits="microMetre" substanceUnits="item"
  volumeUnits="microMetreCubed">
    <notes>
      <body xmlns="http://www.w3.org/1999/xhtml">
        <center>
          <h2>PlazaSur Receptor Dynamics
          </h2>
                </center>
            <p>SBML implementation of the dynamical equation used in PlazaSur.
             <br/>
          Notes on this implementation:
          <ol>
            <li>Assumes that the cell surface area is constant.
            </li>
           <li>The ligand concentration (l) is a constant parameter,
             althought is spatially distributed in the simulations.
           </li>
           <li>The unbound receptor (u) reaches a steady-state value of 52.4.
           </li>
             </ol>
          For further documentation see: http://www.duke.edu/~ccc14/msi-docs/index.html
          <br/>
                </p>
            </body>

    </notes>
    <listOfUnitDefinitions>
      <unitDefinition id="minute">
        <listOfUnits>
          <unit scale="0" exponent="1" multiplier="60" kind="second"/>
        </listOfUnits>
      </unitDefinition>
      <unitDefinition id="microMetre">
        <listOfUnits>
          <unit scale="-6" exponent="1" multiplier="1" kind="metre"/>
        </listOfUnits>
      </unitDefinition>
      <unitDefinition id="microMetreSquared">
        <listOfUnits>
          <unit scale="-6" exponent="2" multiplier="1" kind="metre"/>
        </listOfUnits>
      </unitDefinition>
      <unitDefinition id="microMetreCubed">
        <listOfUnits>
          <unit scale="-6" exponent="3" multiplier="1" kind="metre"/>
        </listOfUnits>
      </unitDefinition>
      <unitDefinition id="perMinute" name="per minute">
        <listOfUnits>
          <unit scale="0" exponent="-1" multiplier="60" kind="second"/>
        </listOfUnits>
      </unitDefinition>
      <unitDefinition id="numPerMicroMetreCubed"
                      name="number per micrometre cubed">
        <listOfUnits>
          <unit scale="0" exponent="1" multiplier="1" kind="item"/>
          <unit scale="-6" exponent="-3" multiplier="1" kind="metre"/>
        </listOfUnits>
      </unitDefinition>
      <unitDefinition id="microMetreCubedPerNum"
                      name="micrometre cubed per number">
        <listOfUnits>
          <unit scale="0" exponent="-1" multiplier="1" kind="item"/>
          <unit scale="-6" exponent="3" multiplier="1" kind="metre"/>
        </listOfUnits>
      </unitDefinition>
      <unitDefinition id="numPerMicroMetreSquaredPerMinute"
                      name="number per micrometre squared per minute">
        <listOfUnits>
          <unit scale="0" exponent="1" multiplier="1" kind="item"/>
          <unit scale="-6" exponent="-2" multiplier="1" kind="metre"/>
          <unit scale="0" exponent="-1" multiplier="60" kind="second"/>
        </listOfUnits>
      </unitDefinition>
    </listOfUnitDefinitions>
    <listOfCompartments>
      <compartment id="cellSurface" constant="true" spatialDimensions="2"
                   units="microMetreSquared" size="314">
        <notes>
          <p xmlns="http://www.w3.org/1999/xhtml">
            The surface area of a 10 micrometer diameter sphere.
          </p>

        </notes>
    </compartment>
    </listOfCompartments>
    <listOfSpecies>
      <species id="u" initialConcentration="100" constant="false"
               hasOnlySubstanceUnits="false" name="unbound receptor"
               boundaryCondition="false" compartment="cellSurface"/>
    </listOfSpecies>
    <listOfParameters>
      <parameter id="a" constant="true" name="gene activity level"
                 value="1" units="dimensionless"/>
      <parameter id="kappa" constant="true" name="affinity"
                 value="0.1" units="microMetreCubedPerNum"/>
      <parameter id="p" constant="true" name="synthesis rate"
                 value="0.1" sboTerm="SBO:0000046"
                 units="numPerMicroMetreSquaredPerMinute"/>
      <parameter id="tau" constant="true"
                 name="internalization rate" value="0.01"
                 units="perMinute"/>
      <parameter id="delta" constant="true" name="degradation rate"
                 value="0.001" sboTerm="SBO:0000356" units="perMinute"/>
      <parameter id="l" constant="true" name="ligand"
                 value="1" units="numPerMicroMetreCubed"/>
    </listOfParameters>
    <listOfReactions>
      <reaction id="synthesis" name="synthesis" reversible="false"
                fast="false" compartment="cellSurface">
        <listOfProducts>
          <speciesReference constant="false" species="u"/>
        </listOfProducts>
        <kineticLaw>
          <math xmlns="http://www.w3.org/1998/Math/MathML">
            <apply>
              <times/>
              <ci> cellSurface
              </ci>
              <apply>
                <times/>
                <ci> p
                </ci>
                <ci> a
                </ci>
              </apply>
            </apply>
          </math>
        </kineticLaw>
      </reaction>
      <reaction id="binding" name="binding and internalization"
                reversible="false" fast="false"
                compartment="cellSurface">
        <listOfReactants>
          <speciesReference constant="false" species="u"/>
        </listOfReactants>
        <kineticLaw>
          <math xmlns="http://www.w3.org/1998/Math/MathML">
            <apply>
              <divide/>
              <apply>
                <times/>
                <apply>
                  <times/>
                  <apply>
                    <times/>
                    <apply>
                      <times/>
                      <ci> cellSurface
                      </ci>
                      <ci> u
                      </ci>
                    </apply>
                    <ci> tau
                    </ci>
                  </apply>
                  <ci> kappa
                  </ci>
                </apply>
                <ci> l
                </ci>
              </apply>
              <apply>
                <plus/>
                <cn type="integer"> 1
                </cn>
                <apply>
                  <times/>
                  <ci> kappa
                  </ci>
                  <ci> l
                  </ci>
                </apply>
              </apply>
            </apply>
          </math>
        </kineticLaw>
      </reaction>
      <reaction id="degradation" name="degradation"
                reversible="false" fast="false"
                compartment="cellSurface">
        <listOfReactants>
          <speciesReference constant="false" species="u"/>
        </listOfReactants>
        <kineticLaw>
          <math xmlns="http://www.w3.org/1998/Math/MathML">
            <apply>
              <times/>
              <apply>
                <times/>
                <ci> delta
                </ci>
                <ci> u
                </ci>
              </apply>
              <ci> cellSurface
              </ci>
            </apply>
          </math>
        </kineticLaw>
      </reaction>
    </listOfReactions>
  </model>
</sbml>