Tutorial: creating a complete model

Table of Contents

• General characteristics of the model to be created
• The program for creating the model
  • Function for creating the model
  • Function for verifying the model
  • Function for writing the model
• Source code

This section presents a complete example program in C++ demonstrating how to create a complete SBML Level 2 model programmatically and then writing the result to a file.

General characteristics of the model to be created

The model presented here is the one described in Section 7.1 of the SBML Level 2 Version 4 specification. It consists of four species identified simply as E, S, ES, and P; one compartment identified as cytosol; and two reactions (one reversible, one irreversible) described by the following equation:

From these equations, it becomes clear that the model also contains three parameters: k_on, k_off, and k_cat.

We will need to assign initial values to these parameters, as well as to the size of the compartment and the initial quantities of the species. Since this model is only meant to illustrate the principles of programming with libSBML, and does not represent an actual biological entity, we hope readers will forgive us for not worrying about justifying the parameter values on the basis of biological plausibility. The table of values for all the quantities we will use is shown below:

Entity in the Model	Initial value	Units
cytosol	1 \(\times\) 10-14	litres
S	1 \(\times\) 10-20	moles
P	0	moles
E	1 \(\times\) 10-21	moles
ES	0	moles
kon	1 000 000	litre/mole/second
koff	0.2	1/second
kcat	0.1	1/second

To add some spice to our boring (but instructive) example, we will use substance units instead of concentration units for the initial values of the species. This is reflected in the use of moles instead of moles/litre as the species quantity units in the table above.

The program for creating the model

To begin, our program has to include certain header files needed by libSBML under C++. For convenience, we will also set the default C++ namespace to std so that we do not have to prefix common C++ functions used throughout our code.

/**
 * File:        example.cpp
 * Description: Create an example model.
 * Author:      Akiya Jouraku
 */
 
#include <iostream>
#include <sbml/SBMLTypes.h>
 
using namespace std;

The next few lines of our program will not be used until the end of the program, but the relevant definitions are best put near the beginning as a principle of good programming. Specifically, this is the definition of a program name and version number, to be inserted by libSBML as XML comments into the final SBML file. We will do something very simple here and define the relevant values as global parameters; other approaches are possible.

//
// These variables are used in writeExampleSBML when writing an SBML
// document.  They are handed to libSBML functions in order to include
// the program information into comments within the SBML file.
//
const static string ProgramName    = "createExample";
const static string ProgramVersion = "1.0.0";

We need a few forward declarations for helper functions. Here they are:

//
// Forward declarations for functions defined below.
//
SBMLDocument* createExampleEnzymaticReaction();
bool validateExampleSBML(SBMLDocument *sbmlDoc);
bool writeExampleSBML(const SBMLDocument *sbmlDoc, const string& filename);

Now we can define a main routine for the program. This essentially acts as scaffolding for calling worker functions and dealing with overall errors.

//
// Main routine
//
int
main (int argc, char *argv[])
{
  SBMLDocument* sbmlDoc = 0;
  bool isValid          = false;
 
  try
  {
    sbmlDoc = createExampleEnzymaticReaction(); 
    isValid = validateExampleSBML(sbmlDoc);
    if (isValid) writeExampleSBML(sbmlDoc, "enzymaticreaction.xml");
 
    // Avoid memory leaks.
    delete sbmlDoc;
 
    // Return non-zero to indicate if a problem occurred.
    if (!isValid) return 1;
  }
  catch (std::bad_alloc& e)
  {
    cerr << e.what() << ": Unable to allocate memory." << endl;
    return 1;
  }
  catch (...)
  {
    cerr << "Unexpected exceptional condition encountered." << endl;
    return 1;
  }
 
  // A 0 return status is the standard Unix/Linux way to say "all ok".
  return 0;
}

The code above has three key functions, createExampleEnzymaticReaction(), validateExampleSBML, and writeExampleSBML(), which do all the real work in the program. Their verbose names hopefully imply their functions. They are described and defined in the next three subsections.

Function for creating the model

We begin this function by creating an SBMLDocument object, creating a Model object within the document object, and setting the Model object's identifier.

/**
 * Create an SBML model represented in "7.1 A Simple example application
 * of SBML" in the SBML Level 2 Version 4 Specification.
 */
SBMLDocument* createExampleEnzymaticReaction()
{
  // Create an SBMLDocument object in Level 2 Version 4 format:
 
  SBMLDocument* sbmlDoc = new SBMLDocument(2, 4);
 
  // Create a Model object inside the SBMLDocument object and set its
  // identifier:
 
  Model* model = sbmlDoc->createModel();
  model->setId("EnzymaticReaction");

Next, we need to define some units in order to properly characterize the quantities used in the model. The libSBML API for unit definitions mirrors the SBML object model, which means there is an object class for UnitDefinition and an object class for the Unit objects that can be placed inside a UnitDefinition object to compose whatever unit is desired. The resulting code for creating units is a little bit tedious because there is a lot to set up, but on the other hand, it is quite straightforward and hopefully easy to follow:

// Temporary pointers (reused more than once below).
 
UnitDefinition* unitdef;
Unit* unit;
 
// (UnitDefinition 1) Create a UnitDefinition object for "per_second".
 
unitdef = model->createUnitDefinition();
unitdef->setId("per_second");
 
// Create an Unit inside the UnitDefinition object above:
 
unit = unitdef->createUnit();
unit->setKind(UNIT_KIND_SECOND);
unit->setExponent(-1);
 
// (UnitDefinition 2) Create a UnitDefinition object for
// "litre_per_mole_per_second".  Note that we can reuse the pointers
// 'unitdef' and 'unit' because the actual UnitDefinition object
// (along with the Unit objects within it) is already attached to
// the Model object. (Refer to the first two lines above.)
  
unitdef = model->createUnitDefinition();
unitdef->setId("litre_per_mole_per_second");
  
// Create individual unit objects that will be put inside
// the UnitDefinition to compose "litre_per_mole_per_second".
 
unit = unitdef->createUnit();
unit->setKind(UNIT_KIND_MOLE);
unit->setExponent(-1);
 
unit = unitdef->createUnit();
unit->setKind(UNIT_KIND_LITRE);
unit->setExponent(1);
 
unit = unitdef->createUnit();
unit->setKind(UNIT_KIND_SECOND);
unit->setExponent(-1);

Before going on, we digress to explain one libSBML API idiom illustrated in the code above and in the rest of the program below. LibSBML provides two main mechanisms for creating objects: class constructors (e.g., Unit() ), and createObject() methods provided by certain object classes such as Model , where Object is another object class such as UnitDefinition . The multiple mechanisms are provided by libSBML for flexibility and to support different use-cases, but they also have different implications for the overall model structure.

In general, the recommended approach is to use the createObject() methods, as in the calls to UnitDefinition::createUnit() above. These methods create an instance and immediately link it into the parent, as well as return a pointer to the object created. By contrast, the other main way of creating an object involves first using the class constructor (e.g., Unit() ) to create an object instance, and then calling an addObject(...) method to add the instance to another object. However, the addObject(...) methods in libSBML copy the object instance given to them. The result is that the caller has to be careful about making changes to the object: changes to the original will not affect the copy added to the model by the addObject(...) method, and moreover, the original must be explicitly deleted at some point or else the program will leak memory. In summary, unless special circumstances apply, most applications are best served using the various createObject() methods.

Back to the example now, we create a compartment object:

// Create a string for the identifier of the compartment.  This
// will be used again later, so it's worth defining a constant here.
 
const string compName = "cytosol";
 
// Create a Compartment object ("cytosol")
 
Compartment* comp = model->createCompartment();
comp->setId(compName);
 
// Set the "size" attribute of the Compartment object.  We are not
// setting the units on the compartment size explicitly, so the units
// of this Compartment object will be the default SBML units of volume,
// which are liters.
 
comp->setSize(1e-14);

Next, we turn to the species definitions in the model. For each species, the program needs to set an identifier, the compartment in which it's located, and optionally the initial quantity of the species in the compartment.

// Temporary pointer (reused more than once below).
 
Species *sp;
 
// Create the Species objects inside the Model object. 
 
// (Species1 ) Create a Species object for species "ES", and set its
// compartment attribute to be the "cytosol" compartment in the model.
 
sp = model->createSpecies();
sp->setCompartment(compName);
sp->setId("ES");
 
// Set the "initialAmount" attribute of the Species object.
//
//  In SBML, the units of a Species object's initial quantity are
//  determined by two attributes, "substanceUnits" and
//  "hasOnlySubstanceUnits", and the "spatialDimensions" attribute
//  of the Compartment object ("cytosol") in which the species
//  object is located.  Here, we are using the default values for
//  "substanceUnits" (which is "mole") and "hasOnlySubstanceUnits"
//  (which is "false").  The compartment in which the species is
//  located uses volume units of liters, so the units of these
//  species (when the species appear in numerical formulas in the
//  model) will be moles/liter.  
//
sp->setInitialAmount(0);
 
// (Species 2) Create a Species object for "P"
 
sp = model->createSpecies();
sp->setCompartment(compName);
sp->setId("P");
sp->setInitialAmount(0);
 
// (Species3) Create a Species object for "S"
 
sp = model->createSpecies();
sp->setCompartment(compName);
sp->setId("S");
sp->setInitialAmount(1e-20);
 
// (Species4) Create a Species object for "E"
 
sp = model->createSpecies();
sp->setCompartment(compName);
sp->setId("E");
sp->setInitialAmount(5e-21 );

It remains only to create objects representing the reactions in the model. In this example, we will create the system as one reversible and one irreversible reaction; alternatively, the system could have been created as three irreversible reactions.

This program has hard-coded versions of the reaction rate formulas only for the purposes of this example. Of course, in a more practical real-life situation, the formulas would not be fixed in this way—one would have methods to read/write/translate expressions originating from different sources (users, files, other programs etc.) to/from the Abstract Syntax Trees (ASTs) used by libSBML.

// Temporary pointers reused below.
 
Reaction* reaction;
SpeciesReference* spr;
Parameter* para;
KineticLaw* kl;
 
// Create Reactant objects inside the Reaction object ("veq"). 
 
reaction = model->createReaction();
reaction->setId("veq");
 
// (Reactant 1) Create a Reactant object that references Species "E"
// in the model.  The object will be created within the reaction in the
// SBML <listOfReactants>.
 
spr = reaction->createReactant();
spr->setSpecies("E");
 
// (Reactant 2) Create a Reactant object that references Species "S"
// in the model.
 
spr = reaction->createReactant();
spr->setSpecies("S");
 
// Create a Product object that references Species "ES" in the model. 
 
spr = reaction->createProduct();
spr->setSpecies("ES");
 
// Create a KineticLaw object inside the Reaction object ("veq").
 
kl  = reaction->createKineticLaw();
 
// Next, we will create an Abstract Syntax Tree representation of the
// following MathML expression:
//
//      <math xmlns=\"http://www.w3.org/1998/Math/MathML\">
//        <apply>
//          <times/> 
//          <ci> cytosol </ci>
//          <apply>
//            <minus/>
//            <apply>
//              <times/>
//              <ci> kon </ci>
//              <ci> E </ci>
//              <ci> S </ci>
//            </apply>
//            <apply>
//              <times/>
//              <ci> koff </ci>
//              <ci> ES </ci>
//            </apply>
//          </apply>
//        </apply>
//      </math>
 
// First, create nodes representing the variables:
 
ASTNode* astCytosol = new ASTNode(AST_NAME);
astCytosol->setName("cytosol");
 
ASTNode* astKon = new ASTNode(AST_NAME);
astKon->setName("kon");
 
ASTNode* astKoff = new ASTNode(AST_NAME);
astKoff->setName("koff");
 
ASTNode* astE = new ASTNode(AST_NAME);
astE->setName("E");
 
ASTNode* astS = new ASTNode(AST_NAME);
astS->setName("S");
 
ASTNode* astES = new ASTNode(AST_NAME);
astES->setName("ES");
 
// Next, we build up the tree bottom-up, beginning with the
// following part:
//
//  <apply>
//    <times/>
//    <ci> koff </ci>
//    <ci> ES </ci>
//  </apply>
 
ASTNode *astTimes1 = new ASTNode(AST_TIMES);
astTimes1->addChild(astKoff);
astTimes1->addChild(astES);
 
// Continuing, we now create nodes representing the following part:
//
//  <apply>
//    <times/>
//    <ci> kon </ci>
//    <ci> E </ci>
//    <ci> S </ci>
//  </apply>
//
// Note: it is worth pointing out that the expression above has 3
// elements inside the <apply>, but libSBML's ASTs are binary trees,
// and only accept two operands.  The following code illustrates how
// to construct such an expression using binary ASTNode objects.
 
ASTNode *astTimes2 = new ASTNode(AST_TIMES);
astTimes2->addChild(astE);
astTimes2->addChild(astS);
 
ASTNode *astTimes = new ASTNode(AST_TIMES);
astTimes->addChild(astKon);
astTimes->addChild(astTimes2);
 
// Now for a node representing the following part:
//
//  <apply>
//    <minus/>
//    <apply>
//      <times/>
//      <ci> kon </ci>
//      <ci> E </ci>
//      <ci> S </ci>
//    </apply>
//    <apply>
//      <times/>
//      <ci> koff </ci>
//      <ci> ES </ci>
//    </apply>
//  </apply>
 
ASTNode *astMinus = new ASTNode(AST_MINUS);
astMinus->addChild(astTimes);
astMinus->addChild(astTimes1);
 
// And finally, the top node, representing the complete expression:
//
// <apply>
//  <times/>
//  <ci> cytosol </ci>
//  <apply>
//    <minus/>
//    <apply>
//      <times/>
//      <ci> kon </ci>
//      <ci> E </ci>
//      <ci> S </ci>
//    </apply>
//    <apply>
//      <times/>
//      <ci> koff </ci>
//      <ci> ES </ci>
//    </apply>
//  </apply>
// </apply>
 
ASTNode* astMath = new ASTNode(AST_TIMES);
astMath->addChild(astCytosol);
astMath->addChild(astMinus);
 
// Now set the <math> element in the KineticLaw object:
 
kl->setMath(astMath);
 
// KineticLaw::setMath(const ASTNode*) sets the math of the
// KineticLaw object to a COPY of the given ASTNode.  The caller
// (i.e., this program) should delete the original ASTNode object to
// avoid a memory leak.
 
delete astMath;
 
// Create local Parameter objects for "kon" and "koff" inside the
// KineticLaw object.  Note how we set the units on the parameter
// values.
 
para = kl->createParameter();
para->setId("kon");
para->setValue(1000000);
para->setUnits("litre_per_mole_per_second");
 
para = kl->createParameter();
para->setId("koff");
para->setValue(0.2);
para->setUnits("per_second");

That was one reaction. One more reaction to go. This time around, for variety and to illustrate another one of libSBML's facilities, we create the AST for the reaction rate formula starting from a string in MathML form.

// (Reaction 2) Create a Reaction object ("vcat") .
 
reaction = model->createReaction();
reaction->setId("vcat");
reaction->setReversible(false);
 
// Create Reactant objects inside the Reaction object ("vcat"). 
 
// (Reactant 1) Create a Reactant object that references Species
// "ES" in the model.
 
spr = reaction->createReactant();
spr->setSpecies("ES");
 
// (Product 1) Create a Product object that references Species "E"
// in the model.
 
spr = reaction->createProduct();
spr->setSpecies("E");
 
// (Product 2) Create a Product object that references Species "P"
// in the model.
 
spr = reaction->createProduct();
spr->setSpecies("P");
 
// Create a KineticLaw object inside the Reaction object ("vcat"). 
 
kl = reaction->createKineticLaw();
 
// To create mathematical expressions, one would typically construct
// an ASTNode tree as the above example code which creates a math of
// another KineticLaw object.  Here, to save some space and
// illustrate another approach of doing it, we will write out the
// formula in MathML form and then use a libSBML convenience
// function to create the ASTNode tree for us.  (However, this is a
// bit dangerous; it's very easy to make mistakes when writing
// MathML by hand, so in a real program, we would not really want to
// do it this way.  We do this here mostly to illustrate this approach.)
 
string mathXMLString = "<math xmlns=\"http://www.w3.org/1998/Math/MathML\">"
                       "  <apply>"
                       "    <times/>"
                       "    <ci> cytosol </ci>"
                       "    <ci> kcat </ci>"
                       "    <ci> ES </ci>"
                       "  </apply>"
                       "</math>";
 
// The following reuses 'astMath' declared earlier in this program.
// Note also how this is careful to delete the structure after handing
// it to setMath().  As noted earlier, setMath() copies the structure
// given to it, so if the caller is finished using it, the structure
// should be deleted in order to avoid leaking memory.
 
astMath = readMathMLFromString(mathXMLString.c_str());
kl->setMath(astMath);
delete astMath;
 
// Create a local Parameter object for "kcat" inside the KineticLaw
// object.
 
para = kl->createParameter();
para->setId("kcat");
para->setValue(0.1);
para->setUnits("per_second");

That concludes all the model creation code. Our createExampleEnzymaticReaction() function now only needs to return the SBMLDocument object that it created.

  // Return the created SBMLDocument object.
  // The returned object must be explicitly deleted by the caller,
  // otherwise a memory leak will happen.
 
  return sbmlDoc;
}

Function for verifying the model

One of the most important features of libSBML is its ability to perform SBML validation to ensure that a model adheres to the SBML specification for whatever Level+Version combination the model uses. Although the example program presented here is short and we can verify its correctness by inspection, in more general programs, it is difficult to ensure that everything in a model is constructed properly. The ability to perform automatic checks then becomes very useful. We therefore take this opportunity to illustrate the use of libSBML's validation and checking facilities.

The rationale behind libSBML's consistency-checking and validation functions is the following. The individual libSBML for creating and setting attributes and other components are generally quite lenient, and allow a caller to compose SBML entities that may not, in the end, represent valid SBML. This allows applications the freedom to do things such as save incomplete models (which is useful when models are being developed over long periods of time), but at the same time, it means that a separate validation step is necessary when a calling program finally wants to finish a complete SBML document.

LibSBML implements verification of SBML in two steps, represented by the two methods SBMLDocument::checkInternalConsistency() and SBMLDocument::checkConsistency(). The former verifies the basic internal consistency and syntax of an SBML document, and the latter implements more elaborate validation rules (both those defined by the SBML specifications, as well as additional rules offered by libSBML). When an application builds up a model programmatically and is finally ready to say "yes, this model is finished", it should call both of these methods to help ensure the correctness and consistency of the finished result.

The situation is different when an application is reading a model from a file or data stream. Upon reading a model, libSBML automatically verifies the syntax and schema consistency of the document. The application then only needs to check the results by interrogating the error log on the SBMLDocument object and calling SBMLDocument::checkConsistency() to check the model against the SBML validation rules.

The following code illustrates a complete sequence of checking a newly-created model using both SBMLDocument::checkInternalConsistency() and SBMLDocument::checkConsistency(), and handling the reporting of problems to the user.

bool validateExampleSBML (SBMLDocument* sbmlDoc)
{
  if ( ! sbmlDoc )
  {
    cerr << "validateExampleSBML: given a null SBML Document" << endl;
    return false;
  }
 
  string consistencyMessages;
  string validationMessages;
  bool noProblems                     = true;
  unsigned int numCheckFailures       = 0;
  unsigned int numConsistencyErrors   = 0;
  unsigned int numConsistencyWarnings = 0;
  unsigned int numValidationErrors    = 0;
  unsigned int numValidationWarnings  = 0;
 
  // LibSBML is lenient when generating models from scratch using the
  // API for creating objects.  Once the whole model is done and before it
  // gets written out, it's important to check that the whole model is in
  // fact complete, consistent and valid.
 
  numCheckFailures = sbmlDoc->checkInternalConsistency();
  if ( numCheckFailures > 0 )
  {
    noProblems = false;
    for ( unsigned int i = 0; i < numCheckFailures; i++ )
    {
      const SBMLError* sbmlErr = sbmlDoc->getError(i);
      if ( sbmlErr->isFatal() || sbmlErr->isError() )
      {
        ++numConsistencyErrors;
      }
      else
      {
        ++numConsistencyWarnings;
      }      
    } 
    ostringstream oss;
    sbmlDoc->printErrors(oss);
    consistencyMessages = oss.str(); 
  }
 
  // If the internal checks fail, it makes little sense to attempt
  // further validation, because the model may be too compromised to
  // be properly interpreted.
 
  if ( numConsistencyErrors > 0 )
  {
    consistencyMessages += "Further validation aborted."; 
  }
  else
  {
    numCheckFailures = sbmlDoc->checkConsistency();
    if ( numCheckFailures > 0 )
    {
      noProblems = false;
      for ( unsigned int i = 0; i < numCheckFailures; i++ )
      {
        const SBMLError* sbmlErr = sbmlDoc->getError(i);
        if ( sbmlErr->isFatal() || sbmlErr->isError() )
        {
          ++numValidationErrors;
        }
        else
        {
          ++numValidationWarnings;
        }      
      } 
      ostringstream oss;
      sbmlDoc->printErrors(oss);
      validationMessages = oss.str(); 
    }
  }
 
  if ( noProblems )
    return true;
  else
  {
    if ( numConsistencyErrors > 0 )
    {
      cout << "ERROR: encountered " << numConsistencyErrors 
           << " consistency error" << (numConsistencyErrors == 1 ? "" : "s")
       << " in model '" << sbmlDoc->getModel()->getId() << "'." << endl;
    }
    if ( numConsistencyWarnings > 0 )
    {
      cout << "Notice: encountered " << numConsistencyWarnings
           << " consistency warning" << (numConsistencyWarnings == 1 ? "" : "s")
       << " in model '" << sbmlDoc->getModel()->getId() << "'." << endl;
    }
    cout << endl << consistencyMessages;
 
    if ( numValidationErrors > 0 )
    {
      cout << "ERROR: encountered " << numValidationErrors
           << " validation error" << (numValidationErrors == 1 ? "" : "s")
       << " in model '" << sbmlDoc->getModel()->getId() << "'." << endl;
    }
    if ( numValidationWarnings > 0 )
    {
      cout << "Notice: encountered " << numValidationWarnings
           << " validation warning" << (numValidationWarnings == 1 ? "" : "s")
       << " in model '" << sbmlDoc->getModel()->getId() << "'." << endl;
    }
    cout << endl << validationMessages;
 
    return ( numConsistencyErrors == 0 && numValidationErrors == 0 );
  }
}

Function for writing the model

The final step of writing out a model is relatively trivial. LibSBML provides a class of object (SBMLWriter) that provides the necessary functionality. All that is left to do is to call it with a desired file name.

/**
 *
 * Writes the given SBMLDocument to the given file.
 *
 */ 
bool writeExampleSBML(const SBMLDocument* sbmlDoc, const string& filename)
{
  SBMLWriter sbmlWriter;
 
  bool result = sbmlWriter.writeSBML(sbmlDoc, filename);
 
  if ( result )
  {
    cout << "Wrote file \"" << filename << "\"" << endl;
    return true;
  }
  else
  {
    cerr << "Failed to write \"" << filename << "\"" << endl;
    return false;
  }
}

Source code

An expanded version of this example is supplied with the libSBML distribution as the file named createExampleSBML.cpp in the directory examples/c++.